Code
library(tidyverse)
library(janitor)
library(readxl)
library(writexl)
library(patchwork)
library(glmmTMB)
library(car)
library(emmeans)
library(nlme)
library(lmerTest)
library(DHARMa)
library(brms)
theme_set(theme_minimal())library(tidyverse)
library(janitor)
library(readxl)
library(writexl)
library(patchwork)
library(glmmTMB)
library(car)
library(emmeans)
library(nlme)
library(lmerTest)
library(DHARMa)
library(brms)
theme_set(theme_minimal())df_combined <- read_csv(here::here("data", "Combined_Master.csv")) %>% clean_names()
df_year_1 <- read_csv(here::here("data", "Year1_Master.csv")) %>% clean_names()
df_year_2 <- read_csv(here::here("data", "Year2_Master.csv")) %>% clean_names()
# Change date format
df_combined <- df_combined %>%
mutate(date = as.Date(date, format="%m_%d_%y")) %>%
select(-starts_with("x"))
df_year_1 <- df_year_1 %>%
mutate(date = as.Date(date, format="%m_%d_%y")) %>%
select(-starts_with("x"))
df_year_2 <- df_year_2 %>%
mutate(date = as.Date(date, format="%m_%d_%y")) %>%
select(-starts_with("x"))
# Save as single Excel File
out <- list("combined" = df_combined, "year_1" = df_year_1, "year_2" = df_year_2)
write_xlsx(out, "data/Master.xlsx")df_combined <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "combined")
df_year_1 <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "year_1")
df_year_2 <- read_xlsx(here::here("data", "Master.xlsx"), sheet = "year_2")df_combined <- df_combined %>%
mutate(
sample_id = as.factor(sample_id),
site = as.factor(site),
block = as.factor(block),
study_year = as.factor(study_year)
)
df_year_1 <- df_year_1 %>%
mutate(
sample_id = as.factor(sample_id),
site = as.factor(site),
block = as.factor(block),
study_year = as.factor(study_year)
)
df_year_2 <- df_year_2 %>%
mutate(
sample_id = as.factor(sample_id),
site = as.factor(site),
block = as.factor(block),
study_year = as.factor(study_year)
)response_vars <- c("nh4", "no3", "mg", "p", "ec", "p_h")
df_combined_LOD_0 <- df_combined %>%
mutate(across(
all_of(response_vars),
~ {
# Replace "LOD" with 0
x <- ifelse(.x == "LOD", 0, .x)
as.numeric(x)
}
))
df_combined_LOD_NA <- df_combined %>%
mutate(across(
all_of(response_vars),
~ {
# Replace "LOD" with NA
x <- ifelse(.x == "LOD", NA, .x)
as.numeric(x)
}
))df_combined %>%
group_by(treatment) %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x == "LOD", na.rm = T)
}
))# A tibble: 6 × 7
treatment nh4 no3 mg p ec p_h
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 MC0 0.774 0.708 0 0.317 0 0
2 MC1 0.843 0.746 0 0.405 0 0
3 MC2 0.848 0.741 0 0.431 0 0
4 PC0 0.850 0.685 0 0.248 0 0
5 PC1 0.776 0.724 0.00971 0.186 0 0
6 PC2 0.836 0.694 0.00962 0.176 0 0
prop_zero_0 <- df_combined_LOD_0 %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x == 0, na.rm = T)
}
)) %>%
pivot_longer(cols = everything(), names_to = "response", values_to = "prop_zero")
prop_zero_NA <- df_combined_LOD_NA %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x == 0, na.rm = T)
}
)) %>%
pivot_longer(cols = everything(), names_to = "response", values_to = "prop_zero")
prop_zero_0 %>%
left_join(prop_zero_NA, by = "response", suffix = c("_0", "_NA")) %>%
arrange(-prop_zero_0)# A tibble: 6 × 3
response prop_zero_0 prop_zero_NA
<chr> <dbl> <dbl>
1 nh4 0.822 0
2 no3 0.716 0
3 p 0.297 0
4 mg 0.00316 0
5 ec 0 0
6 p_h 0 0
prop_zero_0 <- df_combined_LOD_0 %>%
group_by(study_year) %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x == 0, na.rm = T)
}
)) %>%
pivot_longer(cols = -study_year, names_to = "response", values_to = "prop_zero")
prop_zero_NA <- df_combined_LOD_NA %>%
group_by(study_year) %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x == 0, na.rm = T)
}
)) %>%
pivot_longer(cols = -study_year, names_to = "response", values_to = "prop_zero")
prop_zero_0 %>%
left_join(prop_zero_NA, by = c("response", "study_year"), suffix = c("_0", "_NA")) %>%
arrange(study_year, -prop_zero_0)# A tibble: 12 × 4
study_year response prop_zero_0 prop_zero_NA
<fct> <chr> <dbl> <dbl>
1 1 nh4 0.842 0
2 1 no3 0.709 0
3 1 p 0.192 0
4 1 mg 0.0025 0
5 1 ec 0 0
6 1 p_h 0 0
7 2 nh4 0.789 0
8 2 no3 0.728 0
9 2 p 0.484 0
10 2 mg 0.00429 0
11 2 ec 0 0
12 2 p_h 0 0
prop_neg_0 <- df_combined_LOD_0 %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x < 0, na.rm = T)
}
)) %>%
pivot_longer(cols = everything(), names_to = "response", values_to = "prop_negative")
prop_neg_NA <- df_combined_LOD_NA %>%
summarise(across(
all_of(response_vars),
~ {
x <- mean(.x < 0, na.rm = T)
}
)) %>%
pivot_longer(cols = everything(), names_to = "response", values_to = "prop_negative")
prop_neg_0 %>%
left_join(prop_neg_NA, by = "response", suffix = c("_0", "_NA")) %>%
arrange(-prop_negative_0)# A tibble: 6 × 3
response prop_negative_0 prop_negative_NA
<chr> <dbl> <dbl>
1 nh4 0 0
2 no3 0 0
3 mg 0 0
4 p 0 0
5 ec 0 0
6 p_h 0 0
capwords <- function(s, strict = FALSE) {
cap <- function(s) paste(toupper(substring(s, 1, 1)),
{s <- substring(s, 2); if(strict) tolower(s) else s},
sep = "", collapse = " " )
sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))
}
plot_concentration <- function(r, data) {
data %>%
filter(.data[[r]] > 0) %>%
ggplot(aes(x = treatment, y = .data[[r]], fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_point(alpha = 0.3) +
facet_wrap(vars(study_year), scales = "free") +
labs(title = paste0(capwords(r), " Concentration (Non-Zero) by Treatment"), x = "Treatment", y = capwords(r)) +
theme_bw()
}
plots_0 <- map(response_vars, plot_concentration, data = df_combined_LOD_0)
wrap_plots(plots_0, ncol = 2, guides = "collect")# It will be the same since we are plotting non-zeros
# plots_NA <- map(response_vars, plot_concentration, data = df_combined_LOD_NA)
# wrap_plots(plots_NA, ncol = 2, guides = "collect")capwords <- function(s, strict = FALSE) {
cap <- function(s) paste(toupper(substring(s, 1, 1)),
{s <- substring(s, 2); if(strict) tolower(s) else s},
sep = "", collapse = " " )
sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))
}
plot_response <- function(r, data, by) {
data %>%
ggplot(aes(x = date, y = .data[[r]], color = .data[[by]])) +
facet_wrap(vars(study_year), scales = "free") +
geom_point() +
geom_line() +
theme_bw(base_size = 12) +
theme(legend.position = "none") +
labs(
title = paste0(capwords(r), " Over Time"),
x = "Date"
)
}
plots_0 <- map(response_vars, plot_response, data = df_combined_LOD_0, by = "sample_id")
wrap_plots(plots_0, ncol = 2)plots_NA <- map(response_vars, plot_response, data = df_combined_LOD_NA, by = "sample_id")
wrap_plots(plots_NA, ncol = 2)df_model <- df_combined_LOD_0 %>%
filter(study_year == 1) %>% # Can change this to only look at year 2, or comment out to look at both years
mutate(
date = factor(date),
treatment = factor(treatment),
site = factor(site),
block = factor(block),
# Create combined grouping factors
site_block_id = interaction(site, block, drop = TRUE),
trt_site_block_id = interaction(treatment, site, block, drop = TRUE)
)# fit_lme_ec <- lme(
# fixed = ec ~ date * treatment,
# random = list(
# site_block_id = ~ 1, # Random intercept for site:block combinations
# trt_site_block_id = ~ 1 # Random intercept for treatment:site:block combinations
# ),
# data = df_model,
# na.action = na.exclude
# )
#
# Anova(fit_lme_ec)
#
# plot(fitted(fit_lme_ec), residuals(fit_lme_ec),
# xlab = "Fitted Values",
# ylab = "Residuals",
# main = "Residuals vs. Fitted Values")
# abline(h = 0, col = "red", lty = 2)
#
# emmip(fit_lme_ec, treatment ~ date, type = "response", CIs = T) + theme_bw()fit_glmer_gamma_ec <- glmer(
ec ~ date * treatment +
(1 | site_block_id) +
(1 | trt_site_block_id),
data = df_model,
family = Gamma(link = "log"),
control = glmerControl(optimizer = "bobyqa"), # default optimizer did not converge
na.action = na.exclude
)Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
length(par)^2 is not recommended.
summary(fit_glmer_gamma_ec)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Gamma ( log )
Formula: ec ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: df_model
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
147.1 375.2 -16.6 33.1 347
Scaled residuals:
Min 1Q Median 3Q Max
-2.9861 -0.5108 -0.0331 0.4374 4.3962
Random effects:
Groups Name Variance Std.Dev.
trt_site_block_id (Intercept) 0.03329 0.1825
site_block_id (Intercept) 0.02424 0.1557
Residual 0.06125 0.2475
Number of obs: 404, groups: trt_site_block_id, 48; site_block_id, 8
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 0.276677 0.235992 1.172 0.24104
date2024-01-12 -0.029280 0.125662 -0.233 0.81576
date2024-01-23 -0.155530 0.121384 -1.281 0.20009
date2024-02-07 -0.123118 0.121340 -1.015 0.31027
date2024-02-21 -0.077465 0.125634 -0.617 0.53750
date2024-03-06 -0.109187 0.121636 -0.898 0.36937
date2024-03-22 0.018036 0.125707 0.143 0.88592
date2024-04-04 -0.104127 0.131053 -0.795 0.42688
date2024-04-16 -0.268501 0.124017 -2.165 0.03038 *
treatmentMC1 0.351553 0.218209 1.611 0.10716
treatmentMC2 0.345026 0.221048 1.561 0.11856
treatmentPC0 -0.374132 0.218786 -1.710 0.08726 .
treatmentPC1 -0.027530 0.217548 -0.127 0.89930
treatmentPC2 -0.246458 0.229081 -1.076 0.28199
date2024-01-12:treatmentMC1 -0.084005 0.168474 -0.499 0.61804
date2024-01-23:treatmentMC1 0.005998 0.163360 0.037 0.97071
date2024-02-07:treatmentMC1 -0.238434 0.163187 -1.461 0.14399
date2024-02-21:treatmentMC1 -0.364844 0.166539 -2.191 0.02847 *
date2024-03-06:treatmentMC1 -0.469824 0.163802 -2.868 0.00413 **
date2024-03-22:treatmentMC1 -0.450368 0.166540 -2.704 0.00685 **
date2024-04-04:treatmentMC1 -0.471514 0.173754 -2.714 0.00665 **
date2024-04-16:treatmentMC1 -0.609815 0.165547 -3.684 0.00023 ***
date2024-01-12:treatmentMC2 -0.015509 0.169959 -0.091 0.92729
date2024-01-23:treatmentMC2 0.173799 0.166840 1.042 0.29755
date2024-02-07:treatmentMC2 -0.080414 0.167165 -0.481 0.63048
date2024-02-21:treatmentMC2 -0.041221 0.170343 -0.242 0.80879
date2024-03-06:treatmentMC2 0.020059 0.167565 0.120 0.90471
date2024-03-22:treatmentMC2 0.005482 0.170444 0.032 0.97434
date2024-04-04:treatmentMC2 0.006028 0.174387 0.035 0.97243
date2024-04-16:treatmentMC2 -0.147860 0.172094 -0.859 0.39024
date2024-01-12:treatmentPC0 -0.063463 0.166513 -0.381 0.70311
date2024-01-23:treatmentPC0 -0.107943 0.163261 -0.661 0.50851
date2024-02-07:treatmentPC0 0.096256 0.163313 0.589 0.55559
date2024-02-21:treatmentPC0 0.092224 0.166460 0.554 0.57956
date2024-03-06:treatmentPC0 -0.103295 0.163453 -0.632 0.52742
date2024-03-22:treatmentPC0 -0.034863 0.166469 -0.209 0.83411
date2024-04-04:treatmentPC0 -0.161729 0.170559 -0.948 0.34301
date2024-04-16:treatmentPC0 -0.117894 0.165458 -0.713 0.47614
date2024-01-12:treatmentPC1 -0.214503 0.168533 -1.273 0.20310
date2024-01-23:treatmentPC1 -0.182493 0.165601 -1.102 0.27046
date2024-02-07:treatmentPC1 -0.298994 0.163356 -1.830 0.06720 .
date2024-02-21:treatmentPC1 -0.274619 0.166537 -1.649 0.09915 .
date2024-03-06:treatmentPC1 -0.390718 0.163766 -2.386 0.01704 *
date2024-03-22:treatmentPC1 -0.327529 0.166667 -1.965 0.04939 *
date2024-04-04:treatmentPC1 -0.438312 0.170963 -2.564 0.01035 *
date2024-04-16:treatmentPC1 -0.545116 0.170793 -3.192 0.00141 **
date2024-01-12:treatmentPC2 0.051349 0.187448 0.274 0.78413
date2024-01-23:treatmentPC2 -0.108320 0.178830 -0.606 0.54470
date2024-02-07:treatmentPC2 -0.052320 0.179043 -0.292 0.77012
date2024-02-21:treatmentPC2 -0.035076 0.181874 -0.193 0.84707
date2024-03-06:treatmentPC2 -0.190791 0.179275 -1.064 0.28722
date2024-03-22:treatmentPC2 -0.059831 0.181894 -0.329 0.74221
date2024-04-04:treatmentPC2 -0.187168 0.187910 -0.996 0.31922
date2024-04-16:treatmentPC2 -0.265369 0.180586 -1.469 0.14170
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation matrix not shown by default, as p = 54 > 12.
Use print(value, correlation=TRUE) or
vcov(value) if you need it
optimizer (bobyqa) convergence code: 0 (OK)
maxfun < 10 * length(par)^2 is not recommended.
plot(fitted(fit_glmer_gamma_ec), residuals(fit_glmer_gamma_ec),
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)emmip(fit_glmer_gamma_ec, treatment ~ date, type = "response", CIs = T) + theme_bw()fit_lme_ph <- lme(
fixed = p_h ~ as.factor(date) * treatment,
random = list(
site_block_id = ~ 1, # Random intercept for site:block combinations
trt_site_block_id = ~ 1 # Random intercept for treatment:site:block combinations
),
data = df_model,
na.action = na.exclude
)
Anova(fit_lme_ph)Analysis of Deviance Table (Type II tests)
Response: p_h
Chisq Df Pr(>Chisq)
as.factor(date) 59.922 8 4.829e-10 ***
treatment 24.529 5 0.0001717 ***
as.factor(date):treatment 20.744 40 0.9949031
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
S(fit_lme_ph)Linear mixed model fit by REML, Data: df_model
Fixed Effects:
Formula: p_h ~ as.factor(date) * treatment
Estimate Std.Error df t value
(Intercept) 8.0741630 0.1409290 302 57.292
as.factor(date)2024-01-12 0.1075456 0.1441116 302 0.746
as.factor(date)2024-01-23 0.0358370 0.1394633 302 0.257
as.factor(date)2024-02-07 0.0995870 0.1394633 302 0.714
as.factor(date)2024-02-21 -0.1111751 0.1440079 302 -0.772
as.factor(date)2024-03-06 -0.0391630 0.1394633 302 -0.281
as.factor(date)2024-03-22 -0.1610258 0.1441116 302 -1.117
as.factor(date)2024-04-04 -0.0853503 0.1499612 302 -0.569
as.factor(date)2024-04-16 -0.1474028 0.1427589 302 -1.033
treatmentMC1 -0.4743938 0.1612606 35 -2.942
treatmentMC2 -0.4064689 0.1612408 35 -2.521
treatmentPC0 -0.1746424 0.1564647 35 -1.116
treatmentPC1 -0.4088924 0.1564907 35 -2.613
treatmentPC2 -0.3024218 0.1771349 35 -1.707
as.factor(date)2024-01-12:treatmentMC1 0.2313875 0.1975466 302 1.171
as.factor(date)2024-01-23:treatmentMC1 0.4002817 0.1981007 302 2.021
as.factor(date)2024-02-07:treatmentMC1 0.2643938 0.1917050 302 1.379
as.factor(date)2024-02-21:treatmentMC1 0.2601559 0.1950219 302 1.334
as.factor(date)2024-03-06:treatmentMC1 0.2893938 0.1917050 302 1.510
as.factor(date)2024-03-22:treatmentMC1 0.3512565 0.1950973 302 1.800
as.factor(date)2024-04-04:treatmentMC1 0.2825575 0.2025085 302 1.395
as.factor(date)2024-04-16:treatmentMC1 0.2526336 0.1941574 302 1.301
as.factor(date)2024-01-12:treatmentMC2 -0.0064897 0.1950813 302 -0.033
as.factor(date)2024-01-23:treatmentMC2 0.0794538 0.1980666 302 0.401
as.factor(date)2024-02-07:treatmentMC2 0.1427189 0.1916883 302 0.745
as.factor(date)2024-02-21:treatmentMC2 0.0697310 0.1950072 302 0.358
as.factor(date)2024-03-06:treatmentMC2 0.0377189 0.1916883 302 0.197
as.factor(date)2024-03-22:treatmentMC2 0.0395817 0.1950813 302 0.203
as.factor(date)2024-04-04:treatmentMC2 0.1551562 0.1994248 302 0.778
as.factor(date)2024-04-16:treatmentMC2 -0.0322732 0.1971682 302 -0.164
as.factor(date)2024-01-12:treatmentPC0 -0.1445663 0.1911626 302 -0.756
as.factor(date)2024-01-23:treatmentPC0 -0.0991076 0.1876887 302 -0.528
as.factor(date)2024-02-07:treatmentPC0 -0.1791076 0.1876887 302 -0.954
as.factor(date)2024-02-21:treatmentPC0 -0.0870955 0.1910851 302 -0.456
as.factor(date)2024-03-06:treatmentPC0 -0.1316076 0.1876887 302 -0.701
as.factor(date)2024-03-22:treatmentPC0 -0.0959949 0.1911626 302 -0.502
as.factor(date)2024-04-04:treatmentPC0 -0.0754203 0.1956039 302 -0.386
as.factor(date)2024-04-16:treatmentPC0 -0.2146178 0.1901452 302 -1.129
as.factor(date)2024-01-12:treatmentPC1 0.0724544 0.1936173 302 0.374
as.factor(date)2024-01-23:treatmentPC1 0.2155916 0.1901828 302 1.134
as.factor(date)2024-02-07:treatmentPC1 0.2263924 0.1877103 302 1.206
as.factor(date)2024-02-21:treatmentPC1 0.2259045 0.1911024 302 1.182
as.factor(date)2024-03-06:treatmentPC1 0.1713924 0.1877103 302 0.913
as.factor(date)2024-03-22:treatmentPC1 0.3495051 0.1911795 302 1.828
as.factor(date)2024-04-04:treatmentPC1 0.2163297 0.1956153 302 1.106
as.factor(date)2024-04-16:treatmentPC1 0.1436748 0.1964970 302 0.731
as.factor(date)2024-01-12:treatmentPC2 -0.0387264 0.2149142 302 -0.180
as.factor(date)2024-01-23:treatmentPC2 0.1974218 0.2052379 302 0.962
as.factor(date)2024-02-07:treatmentPC2 0.0036718 0.2052379 302 0.018
as.factor(date)2024-02-21:treatmentPC2 0.1006839 0.2083860 302 0.483
as.factor(date)2024-03-06:treatmentPC2 0.0672609 0.2074953 302 0.324
as.factor(date)2024-03-22:treatmentPC2 0.1380346 0.2084665 302 0.662
as.factor(date)2024-04-04:treatmentPC2 0.1142633 0.2147977 302 0.532
as.factor(date)2024-04-16:treatmentPC2 0.0006617 0.2074494 302 0.003
Pr(>|t|)
(Intercept) < 2e-16 ***
as.factor(date)2024-01-12 0.45609
as.factor(date)2024-01-23 0.79738
as.factor(date)2024-02-07 0.47573
as.factor(date)2024-02-21 0.44071
as.factor(date)2024-03-06 0.77905
as.factor(date)2024-03-22 0.26472
as.factor(date)2024-04-04 0.56968
as.factor(date)2024-04-16 0.30265
treatmentMC1 0.00575 **
treatmentMC2 0.01642 *
treatmentPC0 0.27195
treatmentPC1 0.01314 *
treatmentPC2 0.09662 .
as.factor(date)2024-01-12:treatmentMC1 0.24240
as.factor(date)2024-01-23:treatmentMC1 0.04420 *
as.factor(date)2024-02-07:treatmentMC1 0.16886
as.factor(date)2024-02-21:treatmentMC1 0.18321
as.factor(date)2024-03-06:treatmentMC1 0.13220
as.factor(date)2024-03-22:treatmentMC1 0.07279 .
as.factor(date)2024-04-04:treatmentMC1 0.16395
as.factor(date)2024-04-16:treatmentMC1 0.19419
as.factor(date)2024-01-12:treatmentMC2 0.97348
as.factor(date)2024-01-23:treatmentMC2 0.68860
as.factor(date)2024-02-07:treatmentMC2 0.45713
as.factor(date)2024-02-21:treatmentMC2 0.72091
as.factor(date)2024-03-06:treatmentMC2 0.84414
as.factor(date)2024-03-22:treatmentMC2 0.83935
as.factor(date)2024-04-04:treatmentMC2 0.43717
as.factor(date)2024-04-16:treatmentMC2 0.87009
as.factor(date)2024-01-12:treatmentPC0 0.45009
as.factor(date)2024-01-23:treatmentPC0 0.59786
as.factor(date)2024-02-07:treatmentPC0 0.34070
as.factor(date)2024-02-21:treatmentPC0 0.64887
as.factor(date)2024-03-06:treatmentPC0 0.48372
as.factor(date)2024-03-22:treatmentPC0 0.61592
as.factor(date)2024-04-04:treatmentPC0 0.70008
as.factor(date)2024-04-16:treatmentPC0 0.25992
as.factor(date)2024-01-12:treatmentPC1 0.70851
as.factor(date)2024-01-23:treatmentPC1 0.25786
as.factor(date)2024-02-07:treatmentPC1 0.22873
as.factor(date)2024-02-21:treatmentPC1 0.23809
as.factor(date)2024-03-06:treatmentPC1 0.36193
as.factor(date)2024-03-22:treatmentPC1 0.06851 .
as.factor(date)2024-04-04:treatmentPC1 0.26965
as.factor(date)2024-04-16:treatmentPC1 0.46524
as.factor(date)2024-01-12:treatmentPC2 0.85712
as.factor(date)2024-01-23:treatmentPC2 0.33686
as.factor(date)2024-02-07:treatmentPC2 0.98574
as.factor(date)2024-02-21:treatmentPC2 0.62933
as.factor(date)2024-03-06:treatmentPC2 0.74604
as.factor(date)2024-03-22:treatmentPC2 0.50838
as.factor(date)2024-04-04:treatmentPC2 0.59515
as.factor(date)2024-04-16:treatmentPC2 0.99746
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Random effects:
Formula: ~1 | site_block_id
(Intercept)
StdDev: 0.2151
Formula: ~1 | trt_site_block_id %in% site_block_id
(Intercept) Residual
StdDev: 0.1246 0.2419
Number of Observations: 398
Number of Groups:
site_block_id trt_site_block_id %in% site_block_id
8 48
logLik df AIC BIC
-86.80 57 287.60 506.52
plot(fitted(fit_lme_ph), residuals(fit_lme_ph),
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)emmip(fit_lme_ph, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = df_model,
aes(x = date, y = p_h, color = treatment),
alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
theme_bw()Warning: Removed 34 rows containing missing values or values outside the scale range
(`geom_point()`).
response_vars <- c("nh4", "no3", "mg", "p")
lc_model_data <- df_combined %>%
mutate(
nh4_censored = ifelse(nh4 == "LOD", -1, 0),
no3_censored = ifelse(no3 == "LOD", -1, 0),
mg_censored = ifelse(mg == "LOD", -1, 0),
p_censored = ifelse(p == "LOD", -1, 0),
nh4 = ifelse(nh4 == "LOD", nh4_lod, as.numeric(nh4)),
no3 = ifelse(no3 == "LOD", n03_lod, as.numeric(no3)),
mg = ifelse(mg == "LOD", mg_lod, as.numeric(mg)),
p = ifelse(p == "LOD", p_lod, as.numeric(p)),
date = factor(date),
treatment = factor(treatment),
site = factor(site),
block = factor(block),
# Create combined grouping factors
site_block_id = interaction(site, block, drop = TRUE),
trt_site_block_id = interaction(treatment, site, block, drop = TRUE)
)Warning: There were 4 warnings in `mutate()`.
The first warning was:
ℹ In argument: `nh4 = ifelse(nh4 == "LOD", nh4_lod, as.numeric(nh4))`.
Caused by warning in `ifelse()`:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
fit_no3 <- brm(
formula = bf(no3 | cens(no3_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
data = lc_model_data,
family = Gamma(link = "log"),
chains = 4,
warmup = 1000,
iter = 5000,
cores = parallel::detectCores(),
seed = 12345,
file = here::here("models", "fit_no3")
)
fit_nh4 <- brm(
formula = bf(nh4 | cens(nh4_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
data = lc_model_data,
family = Gamma(link = "log"),
chains = 4,
warmup = 1000,
iter = 5000,
cores = parallel::detectCores(),
seed = 12345,
file = here::here("models", "fit_nh4")
)
fit_mg <- brm(
formula = bf(mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
data = lc_model_data,
family = Gamma(link = "log"),
chains = 4,
warmup = 1000,
iter = 5000,
cores = parallel::detectCores(),
seed = 12345,
file = here::here("models", "fit_mg")
)
fit_p <- brm(
formula = bf(p | cens(p_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)),
data = lc_model_data,
family = Gamma(link = "log"),
chains = 4,
warmup = 1000,
iter = 5000,
cores = parallel::detectCores(),
seed = 12345,
file = here::here("models", "fit_p")
)fit_no3 <- read_rds(here::here("models", "fit_no3.rds"))
fit_nh4 <- read_rds(here::here("models", "fit_nh4.rds"))
fit_mg <- read_rds(here::here("models", "fit_mg.rds"))
fit_p <- read_rds(here::here("models", "fit_p.rds"))summary(fit_no3)Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.
Warning: There were 15861 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gamma
Links: mu = log; shape = identity
Formula: no3 | cens(no3_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data (Number of observations: 656)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~site_block_id (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.33 0.62 0.53 3.02 1.03 94 128
~trt_site_block_id (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.41 0.24 1.02 1.98 1.05 67 248
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 6.72 1.26 4.43 9.32 1.62 7
date2024M01M12 -1.56 1.16 -4.26 0.23 1.41 9
date2024M01M23 -3.72 1.45 -6.79 -1.18 1.38 9
date2024M02M07 -6.60 1.41 -9.62 -4.18 1.37 9
date2024M02M21 -194.86 127.24 -502.21 -20.14 1.14 35
date2024M03M06 -178.57 186.66 -679.90 -14.12 2.19 5
date2024M03M22 -5.25 1.49 -8.13 -2.35 1.32 10
date2024M04M04 -254.08 182.26 -644.58 -22.22 1.11 27
date2024M04M16 -167.18 145.07 -580.69 -17.34 1.41 9
date2025M02M08 -4.86 1.30 -7.65 -2.70 1.44 8
date2025M02M15 -5.53 1.26 -8.07 -3.27 1.40 9
date2025M03M04 -6.24 1.39 -9.18 -3.53 1.37 9
date2025M03M16 -7.13 1.37 -10.06 -4.82 1.34 10
date2025M03M27 -346.88 245.23 -705.98 -21.22 1.71 6
date2025M04M11 -236.76 114.45 -422.86 -19.60 1.56 7
treatmentMC1 -1.92 1.50 -5.04 0.68 1.37 9
treatmentMC2 -0.79 1.63 -4.28 2.28 1.44 8
treatmentPC0 -1.74 1.41 -4.51 0.88 1.42 9
treatmentPC1 -0.84 1.55 -4.37 1.84 1.37 9
treatmentPC2 -1.17 1.69 -4.29 2.49 1.43 9
date2024M01M12:treatmentMC1 1.50 1.47 -1.27 4.48 1.19 17
date2024M01M23:treatmentMC1 2.40 1.78 -1.07 5.88 1.22 14
date2024M02M07:treatmentMC1 2.92 1.80 -0.53 6.58 1.20 15
date2024M02M21:treatmentMC1 -117.50 222.02 -552.42 262.24 1.06 67
date2024M03M06:treatmentMC1 -143.62 273.77 -612.54 501.84 1.25 12
date2024M03M22:treatmentMC1 -0.51 1.99 -4.46 3.35 1.17 18
date2024M04M04:treatmentMC1 -94.38 287.42 -585.97 503.19 1.06 77
date2024M04M16:treatmentMC1 -141.80 251.98 -582.65 430.99 1.22 14
date2025M02M08:treatmentMC1 1.47 1.80 -1.85 5.39 1.20 15
date2025M02M15:treatmentMC1 2.54 1.70 -0.77 5.91 1.21 14
date2025M03M04:treatmentMC1 1.98 1.76 -1.43 5.64 1.18 18
date2025M03M16:treatmentMC1 1.34 1.95 -2.25 5.27 1.15 21
date2025M03M27:treatmentMC1 -3.68 311.52 -599.71 583.58 1.36 9
date2025M04M11:treatmentMC1 -122.96 244.67 -558.55 312.67 1.12 23
date2024M01M12:treatmentMC2 0.07 1.57 -2.78 3.47 1.24 13
date2024M01M23:treatmentMC2 -1.75 1.88 -5.18 2.01 1.25 12
date2024M02M07:treatmentMC2 -0.03 1.96 -3.70 3.98 1.22 14
date2024M02M21:treatmentMC2 187.66 127.30 12.93 493.52 1.14 35
date2024M03M06:treatmentMC2 170.50 186.79 5.28 672.38 2.21 5
date2024M03M22:treatmentMC2 -308.11 206.33 -689.71 -16.19 1.02 201
date2024M04M04:treatmentMC2 -76.98 270.95 -558.10 487.48 1.11 25
date2024M04M16:treatmentMC2 161.78 145.04 11.98 575.04 1.42 9
date2025M02M08:treatmentMC2 1.05 1.69 -2.26 4.48 1.28 11
date2025M02M15:treatmentMC2 0.80 1.65 -2.26 4.22 1.25 13
date2025M03M04:treatmentMC2 0.71 1.85 -2.70 4.68 1.23 14
date2025M03M16:treatmentMC2 0.84 1.88 -2.68 4.68 1.21 15
date2025M03M27:treatmentMC2 339.71 245.16 13.96 698.47 1.71 6
date2025M04M11:treatmentMC2 -125.94 223.27 -536.48 286.52 1.11 30
date2024M01M12:treatmentPC0 0.60 1.47 -1.99 3.59 1.22 13
date2024M01M23:treatmentPC0 1.69 1.81 -1.81 5.26 1.28 11
date2024M02M07:treatmentPC0 -306.17 199.86 -678.13 -12.13 1.02 277
date2024M02M21:treatmentPC0 189.64 127.30 14.95 496.61 1.14 35
date2024M03M06:treatmentPC0 -142.05 258.85 -630.23 388.00 1.26 11
date2024M03M22:treatmentPC0 -327.86 207.19 -686.88 -13.66 1.03 100
date2024M04M04:treatmentPC0 -93.27 265.58 -556.34 454.17 1.06 59
date2024M04M16:treatmentPC0 -174.47 262.01 -613.31 395.92 1.24 13
date2025M02M08:treatmentPC0 3.10 1.62 0.34 6.50 1.26 12
date2025M02M15:treatmentPC0 3.71 1.56 0.69 6.83 1.23 14
date2025M03M04:treatmentPC0 3.88 1.81 0.44 7.66 1.22 15
date2025M03M16:treatmentPC0 2.81 1.70 -0.39 6.35 1.22 15
date2025M03M27:treatmentPC0 -0.10 331.60 -587.85 590.65 1.51 7
date2025M04M11:treatmentPC0 232.20 114.61 14.41 417.56 1.57 7
date2024M01M12:treatmentPC1 0.20 1.39 -2.22 3.09 1.28 11
date2024M01M23:treatmentPC1 -1.63 1.85 -5.10 2.03 1.27 11
date2024M02M07:treatmentPC1 0.31 1.82 -2.92 4.01 1.18 16
date2024M02M21:treatmentPC1 187.80 127.19 13.05 494.90 1.14 35
date2024M03M06:treatmentPC1 172.35 186.73 7.62 673.57 2.18 5
date2024M03M22:treatmentPC1 -1.94 2.01 -5.90 2.03 1.15 21
date2024M04M04:treatmentPC1 -61.78 260.22 -531.84 469.68 1.06 57
date2024M04M16:treatmentPC1 160.68 145.03 10.59 573.66 1.41 9
date2025M02M08:treatmentPC1 1.17 1.73 -1.85 4.81 1.23 13
date2025M02M15:treatmentPC1 0.76 1.59 -2.14 3.92 1.22 14
date2025M03M04:treatmentPC1 0.77 1.86 -2.64 4.47 1.18 16
date2025M03M16:treatmentPC1 0.01 1.85 -3.51 3.78 1.14 21
date2025M03M27:treatmentPC1 342.26 245.26 17.22 701.38 1.71 6
date2025M04M11:treatmentPC1 233.82 114.59 16.59 418.52 1.56 7
date2024M01M12:treatmentPC2 0.84 1.78 -3.13 4.49 1.25 13
date2024M01M23:treatmentPC2 0.97 1.95 -3.35 4.88 1.24 14
date2024M02M07:treatmentPC2 -306.38 190.07 -669.51 -14.12 1.03 150
date2024M02M21:treatmentPC2 -155.88 238.12 -584.40 338.19 1.06 65
date2024M03M06:treatmentPC2 173.86 186.78 10.09 675.60 2.17 5
date2024M03M22:treatmentPC2 0.69 2.07 -3.90 4.62 1.24 15
date2024M04M04:treatmentPC2 250.07 182.36 17.75 640.94 1.11 27
date2024M04M16:treatmentPC2 -166.46 245.41 -593.25 357.19 1.19 17
date2025M02M08:treatmentPC2 1.33 1.83 -2.53 4.87 1.35 10
date2025M02M15:treatmentPC2 2.20 1.82 -1.61 5.66 1.38 10
date2025M03M04:treatmentPC2 1.07 2.01 -4.02 4.75 1.29 11
date2025M03M16:treatmentPC2 2.14 2.00 -2.38 6.12 1.30 12
date2025M03M27:treatmentPC2 342.51 245.55 17.33 701.86 1.71 6
date2025M04M11:treatmentPC2 231.90 114.55 14.88 419.06 1.56 7
Tail_ESS
Intercept 22
date2024M01M12 30
date2024M01M23 22
date2024M02M07 28
date2024M02M21 23
date2024M03M06 15
date2024M03M22 105
date2024M04M04 87
date2024M04M16 16
date2025M02M08 20
date2025M02M15 45
date2025M03M04 29
date2025M03M16 41
date2025M03M27 55
date2025M04M11 22
treatmentMC1 38
treatmentMC2 25
treatmentPC0 21
treatmentPC1 28
treatmentPC2 19
date2024M01M12:treatmentMC1 61
date2024M01M23:treatmentMC1 35
date2024M02M07:treatmentMC1 62
date2024M02M21:treatmentMC1 144
date2024M03M06:treatmentMC1 22
date2024M03M22:treatmentMC1 179
date2024M04M04:treatmentMC1 103
date2024M04M16:treatmentMC1 14
date2025M02M08:treatmentMC1 46
date2025M02M15:treatmentMC1 122
date2025M03M04:treatmentMC1 69
date2025M03M16:treatmentMC1 86
date2025M03M27:treatmentMC1 99
date2025M04M11:treatmentMC1 159
date2024M01M12:treatmentMC2 43
date2024M01M23:treatmentMC2 36
date2024M02M07:treatmentMC2 32
date2024M02M21:treatmentMC2 23
date2024M03M06:treatmentMC2 15
date2024M03M22:treatmentMC2 450
date2024M04M04:treatmentMC2 148
date2024M04M16:treatmentMC2 16
date2025M02M08:treatmentMC2 26
date2025M02M15:treatmentMC2 39
date2025M03M04:treatmentMC2 71
date2025M03M16:treatmentMC2 34
date2025M03M27:treatmentMC2 55
date2025M04M11:treatmentMC2 53
date2024M01M12:treatmentPC0 43
date2024M01M23:treatmentPC0 38
date2024M02M07:treatmentPC0 836
date2024M02M21:treatmentPC0 23
date2024M03M06:treatmentPC0 25
date2024M03M22:treatmentPC0 282
date2024M04M04:treatmentPC0 152
date2024M04M16:treatmentPC0 13
date2025M02M08:treatmentPC0 23
date2025M02M15:treatmentPC0 33
date2025M03M04:treatmentPC0 41
date2025M03M16:treatmentPC0 29
date2025M03M27:treatmentPC0 63
date2025M04M11:treatmentPC0 22
date2024M01M12:treatmentPC1 61
date2024M01M23:treatmentPC1 42
date2024M02M07:treatmentPC1 206
date2024M02M21:treatmentPC1 22
date2024M03M06:treatmentPC1 15
date2024M03M22:treatmentPC1 69
date2024M04M04:treatmentPC1 182
date2024M04M16:treatmentPC1 16
date2025M02M08:treatmentPC1 72
date2025M02M15:treatmentPC1 97
date2025M03M04:treatmentPC1 89
date2025M03M16:treatmentPC1 157
date2025M03M27:treatmentPC1 55
date2025M04M11:treatmentPC1 22
date2024M01M12:treatmentPC2 41
date2024M01M23:treatmentPC2 28
date2024M02M07:treatmentPC2 458
date2024M02M21:treatmentPC2 54
date2024M03M06:treatmentPC2 15
date2024M03M22:treatmentPC2 82
date2024M04M04:treatmentPC2 87
date2024M04M16:treatmentPC2 17
date2025M02M08:treatmentPC2 30
date2025M02M15:treatmentPC2 35
date2025M03M04:treatmentPC2 38
date2025M03M16:treatmentPC2 41
date2025M03M27:treatmentPC2 55
date2025M04M11:treatmentPC2 21
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.37 0.04 0.30 0.46 1.01 269 573
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# neff_ratio(fit_no3)
# prior_summary(fit_no3)
# plot(fit, type = "trace", nvariables = 1)
pp_check(fit_no3) + xlim(0, 100)Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 351 rows containing non-finite outside the scale range
(`stat_density()`).
Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_density()`).
pp_check(fit_no3, type = "dens_overlay", ndraws = 100) + xlim(0, 100)Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 3810 rows containing non-finite outside the scale range
(`stat_density()`).
Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_density()`).
# loo(fit_no3)
conditional_effects(fit_no3, effects = "treatment")conditional_effects(fit_no3, effects = "date:treatment")emm_treatments <- emmeans(fit_no3, ~ treatment, type = "response")NOTE: Results may be misleading due to involvement in interactions
emm_treatments treatment response lower.HPD upper.HPD
MC0 0.00e+00 2.22e-16 0.00
MC1 0.00e+00 2.22e-16 0.00
MC2 0.00e+00 2.22e-16 0.00
PC0 0.00e+00 2.22e-16 0.00
PC1 1.41e-08 2.22e-16 0.49
PC2 0.00e+00 2.22e-16 0.00
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
pairs(emm_treatments) contrast ratio lower.HPD upper.HPD
MC0 / MC1 2.41e+18 2.22e-16 2.13e+50
MC0 / MC2 0.00e+00 2.22e-16 7.13e+18
MC0 / PC0 1.22e+19 2.22e-16 1.61e+50
MC0 / PC1 0.00e+00 2.22e-16 0.00e+00
MC0 / PC2 0.00e+00 2.22e-16 1.32e+15
MC1 / MC2 0.00e+00 2.22e-16 0.00e+00
MC1 / PC0 0.00e+00 2.22e-16 3.86e+33
MC1 / PC1 0.00e+00 2.22e-16 0.00e+00
MC1 / PC2 0.00e+00 2.22e-16 0.00e+00
MC2 / PC0 1.10e+27 2.22e-16 2.29e+59
MC2 / PC1 0.00e+00 2.22e-16 1.00e+00
MC2 / PC2 0.00e+00 2.22e-16 4.56e+22
PC0 / PC1 0.00e+00 2.22e-16 0.00e+00
PC0 / PC2 0.00e+00 2.22e-16 1.00e+00
PC1 / PC2 8.26e+18 4.07e-16 2.71e+38
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
emmip(fit_no3, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data,
aes(x = date, y = no3, color = as.factor(no3_censored)),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
# ylim(0, 1000) +
theme_bw()Warning: Removed 64 rows containing missing values or values outside the scale range
(`geom_point()`).
summary(fit_nh4)Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.
Warning: There were 15995 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gamma
Links: mu = log; shape = identity
Formula: nh4 | cens(nh4_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data (Number of observations: 663)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~site_block_id (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.51 0.41 0.01 1.56 1.05 47 88
~trt_site_block_id (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.86 0.36 1.20 2.74 1.34 9 21
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 0.27 0.98 -1.45 2.35 1.23 15
date2024M01M12 -1.64 1.15 -4.14 0.37 1.11 30
date2024M01M23 -110.65 142.37 -410.51 -6.52 1.68 6
date2024M02M07 -246.28 182.22 -671.07 -15.53 1.49 8
date2024M02M21 -449.79 193.63 -696.20 -48.32 1.45 8
date2024M03M06 -0.92 1.32 -3.86 1.42 1.09 38
date2024M03M22 -3.46 1.66 -6.67 -0.15 1.07 54
date2024M04M04 -0.96 1.00 -3.11 0.85 1.05 57
date2024M04M16 -3.16 1.14 -5.42 -0.94 1.08 45
date2025M02M08 -1.89 1.07 -3.89 0.31 1.22 16
date2025M02M15 -3.30 1.49 -6.33 -0.65 1.22 14
date2025M03M04 -1.91 1.49 -4.39 1.92 1.22 15
date2025M03M16 -3.03 1.18 -5.51 -0.84 1.12 25
date2025M03M27 -1.25 1.08 -3.28 1.02 1.13 34
date2025M04M11 -47.63 61.30 -183.79 -3.44 1.87 6
treatmentMC1 0.97 1.55 -1.71 4.46 1.18 16
treatmentMC2 -0.41 1.30 -3.17 1.87 1.15 23
treatmentPC0 0.59 1.22 -1.97 2.96 1.09 45
treatmentPC1 0.28 1.52 -2.88 3.12 1.24 16
treatmentPC2 -1.09 1.68 -4.40 2.27 1.17 19
date2024M01M12:treatmentMC1 -1.64 1.80 -5.12 2.10 1.15 21
date2024M01M23:treatmentMC1 106.70 142.45 2.75 405.97 1.69 6
date2024M02M07:treatmentMC1 -82.83 252.90 -563.34 404.13 1.24 13
date2024M02M21:treatmentMC1 149.14 266.31 -424.06 598.51 1.33 10
date2024M03M06:treatmentMC1 -414.26 193.97 -699.10 -32.68 1.25 12
date2024M03M22:treatmentMC1 -229.75 196.40 -672.44 -9.15 1.12 27
date2024M04M04:treatmentMC1 -2.13 1.69 -5.69 0.96 1.07 56
date2024M04M16:treatmentMC1 -300.57 197.17 -681.31 -14.27 1.01 168
date2025M02M08:treatmentMC1 -1.41 1.61 -4.79 1.51 1.09 34
date2025M02M15:treatmentMC1 -327.04 199.96 -680.54 -14.20 1.04 88
date2025M03M04:treatmentMC1 -1.13 2.10 -5.90 2.40 1.32 10
date2025M03M16:treatmentMC1 -338.93 205.82 -685.30 -16.92 1.13 30
date2025M03M27:treatmentMC1 -335.21 203.80 -673.01 -15.45 1.12 25
date2025M04M11:treatmentMC1 -302.47 220.82 -671.75 116.09 1.14 21
date2024M01M12:treatmentMC2 0.49 1.71 -2.93 3.98 1.15 31
date2024M01M23:treatmentMC2 107.64 142.51 3.97 408.47 1.68 6
date2024M02M07:treatmentMC2 -70.20 300.91 -601.10 545.05 1.46 8
date2024M02M21:treatmentMC2 107.52 279.08 -453.36 596.20 1.26 11
date2024M03M06:treatmentMC2 -355.19 210.25 -687.27 -16.92 1.12 23
date2024M03M22:treatmentMC2 -361.38 206.46 -685.03 -19.29 1.02 118
date2024M04M04:treatmentMC2 -1.02 1.67 -4.43 2.16 1.09 45
date2024M04M16:treatmentMC2 -343.27 199.20 -680.85 -18.64 1.05 89
date2025M02M08:treatmentMC2 0.23 1.46 -2.56 3.37 1.17 17
date2025M02M15:treatmentMC2 -315.23 189.98 -676.07 -17.03 1.02 147
date2025M03M04:treatmentMC2 -0.09 2.06 -4.96 3.70 1.19 18
date2025M03M16:treatmentMC2 -330.48 198.16 -684.94 -13.04 1.01 139
date2025M03M27:treatmentMC2 0.37 1.50 -2.47 3.26 1.14 25
date2025M04M11:treatmentMC2 46.89 61.36 1.55 182.38 1.85 6
date2024M01M12:treatmentPC0 -340.87 203.12 -690.13 -21.12 1.10 31
date2024M01M23:treatmentPC0 -214.39 251.67 -641.08 320.56 1.26 12
date2024M02M07:treatmentPC0 -94.06 268.29 -572.14 456.48 1.22 13
date2024M02M21:treatmentPC0 115.00 263.95 -383.78 622.94 1.18 16
date2024M03M06:treatmentPC0 -347.83 196.30 -680.71 -20.69 1.02 195
date2024M03M22:treatmentPC0 -317.88 221.06 -696.89 -11.40 1.09 36
date2024M04M04:treatmentPC0 -3.11 1.69 -6.55 0.26 1.04 84
date2024M04M16:treatmentPC0 -0.64 1.79 -3.98 2.95 1.11 42
date2025M02M08:treatmentPC0 0.11 1.70 -2.74 3.59 1.20 16
date2025M02M15:treatmentPC0 -0.35 2.15 -4.26 4.13 1.26 13
date2025M03M04:treatmentPC0 -0.35 2.07 -5.13 3.50 1.23 13
date2025M03M16:treatmentPC0 -1.21 1.85 -4.48 2.81 1.08 47
date2025M03M27:treatmentPC0 -0.79 1.72 -4.07 2.67 1.11 33
date2025M04M11:treatmentPC0 46.03 61.11 1.16 181.12 1.85 6
date2024M01M12:treatmentPC1 -2.28 1.79 -5.52 1.72 1.07 50
date2024M01M23:treatmentPC1 108.06 142.48 4.36 408.26 1.69 6
date2024M02M07:treatmentPC1 -67.53 276.98 -558.40 575.11 1.32 11
date2024M02M21:treatmentPC1 445.63 193.65 44.10 692.41 1.45 8
date2024M03M06:treatmentPC1 -330.59 194.33 -674.70 -17.14 1.01 68
date2024M03M22:treatmentPC1 -300.04 195.33 -677.19 -17.34 1.06 103
date2024M04M04:treatmentPC1 -0.99 1.77 -4.44 2.52 1.05 64
date2024M04M16:treatmentPC1 1.52 1.85 -2.07 4.99 1.08 49
date2025M02M08:treatmentPC1 -0.49 1.86 -3.90 3.09 1.13 30
date2025M02M15:treatmentPC1 0.26 2.01 -3.25 4.40 1.20 16
date2025M03M04:treatmentPC1 0.93 2.04 -3.94 4.70 1.19 18
date2025M03M16:treatmentPC1 0.96 1.89 -2.70 4.86 1.10 30
date2025M03M27:treatmentPC1 0.90 1.80 -2.56 4.34 1.10 36
date2025M04M11:treatmentPC1 48.23 61.39 3.57 184.31 1.92 6
date2024M01M12:treatmentPC2 -1.58 1.85 -5.50 1.81 1.12 24
date2024M01M23:treatmentPC2 109.03 142.38 3.50 409.48 1.68 6
date2024M02M07:treatmentPC2 243.45 182.27 13.00 668.96 1.49 8
date2024M02M21:treatmentPC2 76.69 271.22 -564.20 547.48 1.23 13
date2024M03M06:treatmentPC2 -348.29 204.86 -685.51 -17.95 1.02 103
date2024M03M22:treatmentPC2 -278.55 190.88 -674.71 -14.59 1.10 38
date2024M04M04:treatmentPC2 -1.89 2.03 -5.83 2.15 1.08 54
date2024M04M16:treatmentPC2 -314.52 188.60 -668.68 -16.19 1.03 153
date2025M02M08:treatmentPC2 -0.28 1.92 -4.01 3.25 1.21 17
date2025M02M15:treatmentPC2 -341.83 200.21 -684.05 -13.71 1.03 97
date2025M03M04:treatmentPC2 -0.81 2.10 -5.94 2.73 1.13 26
date2025M03M16:treatmentPC2 2.10 1.84 -1.91 5.68 1.23 13
date2025M03M27:treatmentPC2 0.60 1.89 -3.14 4.29 1.18 19
date2025M04M11:treatmentPC2 48.59 61.07 3.87 184.03 1.87 6
Tail_ESS
Intercept 79
date2024M01M12 120
date2024M01M23 21
date2024M02M07 15
date2024M02M21 22
date2024M03M06 125
date2024M03M22 100
date2024M04M04 96
date2024M04M16 88
date2025M02M08 52
date2025M02M15 67
date2025M03M04 15
date2025M03M16 118
date2025M03M27 120
date2025M04M11 20
treatmentMC1 50
treatmentMC2 44
treatmentPC0 88
treatmentPC1 71
treatmentPC2 36
date2024M01M12:treatmentMC1 134
date2024M01M23:treatmentMC1 20
date2024M02M07:treatmentMC1 65
date2024M02M21:treatmentMC1 71
date2024M03M06:treatmentMC1 98
date2024M03M22:treatmentMC1 164
date2024M04M04:treatmentMC1 105
date2024M04M16:treatmentMC1 387
date2025M02M08:treatmentMC1 97
date2025M02M15:treatmentMC1 358
date2025M03M04:treatmentMC1 18
date2025M03M16:treatmentMC1 229
date2025M03M27:treatmentMC1 321
date2025M04M11:treatmentMC1 36
date2024M01M12:treatmentMC2 63
date2024M01M23:treatmentMC2 20
date2024M02M07:treatmentMC2 14
date2024M02M21:treatmentMC2 47
date2024M03M06:treatmentMC2 296
date2024M03M22:treatmentMC2 361
date2024M04M04:treatmentMC2 118
date2024M04M16:treatmentMC2 461
date2025M02M08:treatmentMC2 89
date2025M02M15:treatmentMC2 356
date2025M03M04:treatmentMC2 41
date2025M03M16:treatmentMC2 342
date2025M03M27:treatmentMC2 124
date2025M04M11:treatmentMC2 19
date2024M01M12:treatmentPC0 178
date2024M01M23:treatmentPC0 31
date2024M02M07:treatmentPC0 24
date2024M02M21:treatmentPC0 58
date2024M03M06:treatmentPC0 480
date2024M03M22:treatmentPC0 227
date2024M04M04:treatmentPC0 137
date2024M04M16:treatmentPC0 96
date2025M02M08:treatmentPC0 108
date2025M02M15:treatmentPC0 34
date2025M03M04:treatmentPC0 13
date2025M03M16:treatmentPC0 122
date2025M03M27:treatmentPC0 118
date2025M04M11:treatmentPC0 20
date2024M01M12:treatmentPC1 118
date2024M01M23:treatmentPC1 21
date2024M02M07:treatmentPC1 12
date2024M02M21:treatmentPC1 22
date2024M03M06:treatmentPC1 247
date2024M03M22:treatmentPC1 300
date2024M04M04:treatmentPC1 138
date2024M04M16:treatmentPC1 69
date2025M02M08:treatmentPC1 125
date2025M02M15:treatmentPC1 67
date2025M03M04:treatmentPC1 20
date2025M03M16:treatmentPC1 130
date2025M03M27:treatmentPC1 46
date2025M04M11:treatmentPC1 18
date2024M01M12:treatmentPC2 107
date2024M01M23:treatmentPC2 20
date2024M02M07:treatmentPC2 15
date2024M02M21:treatmentPC2 40
date2024M03M06:treatmentPC2 489
date2024M03M22:treatmentPC2 477
date2024M04M04:treatmentPC2 120
date2024M04M16:treatmentPC2 390
date2025M02M08:treatmentPC2 52
date2025M02M15:treatmentPC2 400
date2025M03M04:treatmentPC2 29
date2025M03M16:treatmentPC2 30
date2025M03M27:treatmentPC2 73
date2025M04M11:treatmentPC2 18
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.30 0.05 0.21 0.40 1.08 47 212
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# neff_ratio(fit_nh4)
# prior_summary(fit_nh4)
# plot(fit, type = "trace", nvariables = 1)
pp_check(fit_nh4) + xlim(0, 10)Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 62 rows containing non-finite outside the scale range
(`stat_density()`).
pp_check(fit_nh4, type = "dens_overlay", ndraws = 100) + xlim(0, 10)Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 726 rows containing non-finite outside the scale range
(`stat_density()`).
# loo(fit_nh4)
conditional_effects(fit_nh4, effects = "treatment")conditional_effects(fit_nh4, effects = "date:treatment")emm_treatments <- emmeans(fit_nh4, ~ treatment, type = "response")NOTE: Results may be misleading due to involvement in interactions
emm_treatments treatment response lower.HPD upper.HPD
MC0 2.2204e-16 2.2204e-16 2.1125e-13
MC1 2.2204e-16 2.2204e-16 2.2200e-16
MC2 2.2204e-16 2.2204e-16 2.2200e-16
PC0 2.2204e-16 2.2204e-16 2.2200e-16
PC1 2.2204e-16 2.2204e-16 7.5358e-13
PC2 2.2204e-16 2.2204e-16 2.2200e-16
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
pairs(emm_treatments) contrast ratio lower.HPD upper.HPD
MC0 / MC1 3.18e+60 148877 1.60e+86
MC0 / MC2 1.10e+43 0 2.11e+76
MC0 / PC0 1.74e+34 0 8.32e+56
MC0 / PC1 6.70e+01 0 2.13e+25
MC0 / PC2 5.49e+23 0 3.70e+48
MC1 / MC2 0.00e+00 0 2.32e+22
MC1 / PC0 0.00e+00 0 1.82e+04
MC1 / PC1 0.00e+00 0 0.00e+00
MC1 / PC2 0.00e+00 0 0.00e+00
MC2 / PC0 0.00e+00 0 1.66e+30
MC2 / PC1 0.00e+00 0 0.00e+00
MC2 / PC2 0.00e+00 0 1.14e+14
PC0 / PC1 0.00e+00 0 0.00e+00
PC0 / PC2 0.00e+00 0 5.06e+24
PC1 / PC2 5.85e+21 0 9.11e+47
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
plot(emm_treatments) + theme_bw()emmip(fit_nh4, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data,
aes(x = date, y = nh4, color = as.factor(mg_censored)),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
# ylim(0, 10) +
theme_bw()Warning: Removed 57 rows containing missing values or values outside the scale range
(`geom_point()`).
summary(fit_mg)Warning: There were 149 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gamma
Links: mu = log; shape = identity
Formula: mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data (Number of observations: 633)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~site_block_id (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.40 0.18 0.17 0.87 1.00 2032 743
~trt_site_block_id (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.38 0.05 0.29 0.50 1.00 6068 9749
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 3.90 0.30 3.29 4.45 1.00 869
date2024M01M12 0.20 0.24 -0.28 0.68 1.00 2141
date2024M01M23 0.25 0.24 -0.23 0.71 1.00 1933
date2024M02M07 -0.01 0.24 -0.48 0.46 1.00 2122
date2024M02M21 0.12 0.25 -0.36 0.62 1.00 2153
date2024M03M06 0.19 0.25 -0.30 0.66 1.00 1806
date2024M03M22 0.38 0.24 -0.09 0.86 1.00 1848
date2024M04M04 0.48 0.26 -0.02 0.98 1.00 2203
date2024M04M16 0.34 0.24 -0.14 0.81 1.00 2128
date2025M02M08 -0.66 0.24 -1.14 -0.19 1.00 1980
date2025M02M15 -0.48 0.25 -0.97 0.00 1.00 2085
date2025M03M04 -0.35 0.26 -0.85 0.16 1.00 2299
date2025M03M16 -0.14 0.24 -0.61 0.33 1.00 2104
date2025M03M27 0.03 0.26 -0.47 0.53 1.00 2299
date2025M04M11 0.21 0.27 -0.31 0.75 1.00 1922
treatmentMC1 0.86 0.31 0.25 1.47 1.00 2513
treatmentMC2 0.58 0.32 -0.04 1.22 1.00 2455
treatmentPC0 0.01 0.32 -0.62 0.62 1.00 2577
treatmentPC1 0.48 0.31 -0.12 1.09 1.00 2187
treatmentPC2 0.10 0.35 -0.57 0.79 1.00 1728
date2024M01M12:treatmentMC1 -0.27 0.33 -0.92 0.38 1.00 2988
date2024M01M23:treatmentMC1 -0.32 0.33 -0.96 0.33 1.00 2830
date2024M02M07:treatmentMC1 -0.55 0.32 -1.20 0.07 1.00 2961
date2024M02M21:treatmentMC1 -0.73 0.33 -1.39 -0.09 1.00 2731
date2024M03M06:treatmentMC1 -0.81 0.33 -1.47 -0.15 1.00 2874
date2024M03M22:treatmentMC1 -0.80 0.33 -1.44 -0.15 1.00 2825
date2024M04M04:treatmentMC1 -0.92 0.35 -1.60 -0.23 1.00 2874
date2024M04M16:treatmentMC1 -1.01 0.33 -1.66 -0.37 1.00 2227
date2025M02M08:treatmentMC1 0.01 0.34 -0.67 0.66 1.00 2500
date2025M02M15:treatmentMC1 -0.07 0.33 -0.72 0.57 1.00 2770
date2025M03M04:treatmentMC1 -0.10 0.34 -0.78 0.56 1.00 3144
date2025M03M16:treatmentMC1 -0.47 0.32 -1.12 0.16 1.00 2787
date2025M03M27:treatmentMC1 -0.50 0.34 -1.19 0.16 1.00 2750
date2025M04M11:treatmentMC1 -0.40 0.36 -1.11 0.29 1.00 3045
date2024M01M12:treatmentMC2 -0.04 0.33 -0.69 0.62 1.00 2875
date2024M01M23:treatmentMC2 0.06 0.33 -0.59 0.70 1.00 2775
date2024M02M07:treatmentMC2 -0.19 0.33 -0.84 0.45 1.00 2798
date2024M02M21:treatmentMC2 -0.14 0.34 -0.82 0.51 1.00 2710
date2024M03M06:treatmentMC2 -0.10 0.34 -0.76 0.57 1.00 2825
date2024M03M22:treatmentMC2 -0.11 0.33 -0.77 0.53 1.00 2511
date2024M04M04:treatmentMC2 -0.19 0.34 -0.86 0.48 1.00 2886
date2024M04M16:treatmentMC2 -0.37 0.34 -1.04 0.29 1.00 2516
date2025M02M08:treatmentMC2 0.19 0.35 -0.49 0.86 1.00 2463
date2025M02M15:treatmentMC2 0.66 0.34 -0.01 1.31 1.00 2887
date2025M03M04:treatmentMC2 0.55 0.35 -0.15 1.23 1.00 3094
date2025M03M16:treatmentMC2 0.21 0.33 -0.44 0.85 1.00 2700
date2025M03M27:treatmentMC2 0.12 0.34 -0.56 0.79 1.00 3062
date2025M04M11:treatmentMC2 0.12 0.39 -0.64 0.88 1.00 3304
date2024M01M12:treatmentPC0 -0.31 0.33 -0.95 0.35 1.00 2992
date2024M01M23:treatmentPC0 -0.56 0.32 -1.19 0.07 1.00 2745
date2024M02M07:treatmentPC0 -0.40 0.33 -1.05 0.23 1.00 2922
date2024M02M21:treatmentPC0 -0.37 0.33 -1.02 0.28 1.00 3116
date2024M03M06:treatmentPC0 -0.36 0.33 -1.00 0.30 1.00 2483
date2024M03M22:treatmentPC0 -0.54 0.33 -1.19 0.11 1.00 2706
date2024M04M04:treatmentPC0 -0.59 0.34 -1.25 0.09 1.00 2953
date2024M04M16:treatmentPC0 -0.51 0.33 -1.17 0.13 1.00 3257
date2025M02M08:treatmentPC0 -0.59 0.39 -1.35 0.17 1.00 3075
date2025M02M15:treatmentPC0 -0.55 0.33 -1.21 0.10 1.00 2969
date2025M03M04:treatmentPC0 -0.46 0.36 -1.17 0.26 1.00 3014
date2025M03M16:treatmentPC0 -0.81 0.34 -1.47 -0.14 1.00 3140
date2025M03M27:treatmentPC0 -0.67 0.38 -1.39 0.08 1.00 3335
date2025M04M11:treatmentPC0 -0.53 0.37 -1.25 0.19 1.00 2247
date2024M01M12:treatmentPC1 -0.46 0.33 -1.11 0.19 1.00 2995
date2024M01M23:treatmentPC1 -0.60 0.33 -1.23 0.05 1.00 2892
date2024M02M07:treatmentPC1 -0.91 0.32 -1.54 -0.29 1.00 2752
date2024M02M21:treatmentPC1 -0.84 0.33 -1.48 -0.21 1.00 2476
date2024M03M06:treatmentPC1 -0.73 0.33 -1.37 -0.08 1.00 2880
date2024M03M22:treatmentPC1 -0.82 0.33 -1.46 -0.18 1.00 2772
date2024M04M04:treatmentPC1 -0.95 0.34 -1.61 -0.28 1.00 2979
date2024M04M16:treatmentPC1 -1.03 0.34 -1.70 -0.37 1.00 2945
date2025M02M08:treatmentPC1 -1.01 0.33 -1.66 -0.36 1.00 2785
date2025M02M15:treatmentPC1 -1.05 0.34 -1.72 -0.39 1.00 2882
date2025M03M04:treatmentPC1 -1.12 0.34 -1.80 -0.45 1.00 3182
date2025M03M16:treatmentPC1 -1.05 0.34 -1.71 -0.39 1.00 3087
date2025M03M27:treatmentPC1 -1.04 0.35 -1.72 -0.34 1.00 3212
date2025M04M11:treatmentPC1 -1.23 0.40 -2.00 -0.44 1.00 3125
date2024M01M12:treatmentPC2 -0.29 0.36 -1.00 0.43 1.00 2778
date2024M01M23:treatmentPC2 -0.55 0.36 -1.25 0.14 1.00 2333
date2024M02M07:treatmentPC2 -0.78 0.36 -1.47 -0.07 1.00 2159
date2024M02M21:treatmentPC2 -0.61 0.36 -1.31 0.11 1.00 2491
date2024M03M06:treatmentPC2 -0.58 0.37 -1.30 0.13 1.00 1856
date2024M03M22:treatmentPC2 -0.55 0.36 -1.26 0.15 1.00 2027
date2024M04M04:treatmentPC2 -0.62 0.37 -1.34 0.12 1.00 2706
date2024M04M16:treatmentPC2 -0.68 0.36 -1.39 0.02 1.00 3141
date2025M02M08:treatmentPC2 -0.65 0.37 -1.37 0.07 1.00 2916
date2025M02M15:treatmentPC2 -0.36 0.36 -1.07 0.35 1.00 2638
date2025M03M04:treatmentPC2 -0.22 0.38 -0.96 0.52 1.00 2541
date2025M03M16:treatmentPC2 -0.40 0.36 -1.11 0.30 1.00 2830
date2025M03M27:treatmentPC2 -0.38 0.38 -1.14 0.36 1.00 2334
date2025M04M11:treatmentPC2 -0.56 0.43 -1.38 0.30 1.00 2570
Tail_ESS
Intercept 489
date2024M01M12 4611
date2024M01M23 4161
date2024M02M07 4299
date2024M02M21 5113
date2024M03M06 5267
date2024M03M22 5165
date2024M04M04 5257
date2024M04M16 5252
date2025M02M08 4562
date2025M02M15 5271
date2025M03M04 5281
date2025M03M16 5087
date2025M03M27 5755
date2025M04M11 5094
treatmentMC1 5223
treatmentMC2 6142
treatmentPC0 6150
treatmentPC1 5902
treatmentPC2 1876
date2024M01M12:treatmentMC1 7291
date2024M01M23:treatmentMC1 6716
date2024M02M07:treatmentMC1 6611
date2024M02M21:treatmentMC1 5988
date2024M03M06:treatmentMC1 6631
date2024M03M22:treatmentMC1 6946
date2024M04M04:treatmentMC1 6233
date2024M04M16:treatmentMC1 5203
date2025M02M08:treatmentMC1 6029
date2025M02M15:treatmentMC1 5923
date2025M03M04:treatmentMC1 6228
date2025M03M16:treatmentMC1 5853
date2025M03M27:treatmentMC1 7492
date2025M04M11:treatmentMC1 7125
date2024M01M12:treatmentMC2 7081
date2024M01M23:treatmentMC2 6424
date2024M02M07:treatmentMC2 6043
date2024M02M21:treatmentMC2 6387
date2024M03M06:treatmentMC2 7069
date2024M03M22:treatmentMC2 5811
date2024M04M04:treatmentMC2 6399
date2024M04M16:treatmentMC2 7370
date2025M02M08:treatmentMC2 6477
date2025M02M15:treatmentMC2 6186
date2025M03M04:treatmentMC2 6360
date2025M03M16:treatmentMC2 6414
date2025M03M27:treatmentMC2 7132
date2025M04M11:treatmentMC2 6804
date2024M01M12:treatmentPC0 6866
date2024M01M23:treatmentPC0 7129
date2024M02M07:treatmentPC0 6789
date2024M02M21:treatmentPC0 6619
date2024M03M06:treatmentPC0 5943
date2024M03M22:treatmentPC0 7353
date2024M04M04:treatmentPC0 6545
date2024M04M16:treatmentPC0 7638
date2025M02M08:treatmentPC0 8050
date2025M02M15:treatmentPC0 6781
date2025M03M04:treatmentPC0 7508
date2025M03M16:treatmentPC0 6850
date2025M03M27:treatmentPC0 7941
date2025M04M11:treatmentPC0 5367
date2024M01M12:treatmentPC1 6645
date2024M01M23:treatmentPC1 6957
date2024M02M07:treatmentPC1 7158
date2024M02M21:treatmentPC1 6270
date2024M03M06:treatmentPC1 7322
date2024M03M22:treatmentPC1 6564
date2024M04M04:treatmentPC1 6755
date2024M04M16:treatmentPC1 7499
date2025M02M08:treatmentPC1 6374
date2025M02M15:treatmentPC1 7465
date2025M03M04:treatmentPC1 7572
date2025M03M16:treatmentPC1 7202
date2025M03M27:treatmentPC1 7874
date2025M04M11:treatmentPC1 8764
date2024M01M12:treatmentPC2 6428
date2024M01M23:treatmentPC2 6396
date2024M02M07:treatmentPC2 6143
date2024M02M21:treatmentPC2 6478
date2024M03M06:treatmentPC2 3113
date2024M03M22:treatmentPC2 4921
date2024M04M04:treatmentPC2 7305
date2024M04M16:treatmentPC2 6883
date2025M02M08:treatmentPC2 7278
date2025M02M15:treatmentPC2 7244
date2025M03M04:treatmentPC2 6715
date2025M03M16:treatmentPC2 6105
date2025M03M27:treatmentPC2 3233
date2025M04M11:treatmentPC2 7080
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 5.78 0.35 5.11 6.48 1.00 13355 11217
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# neff_ratio(fit_mg)
# prior_summary(fit_mg)
# plot(fit, type = "trace")
pp_check(fit_mg)Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pp_check(fit_mg, type = "dens_overlay", ndraws = 100)Warning: Censored responses are not shown in 'pp_check'.
# loo(fit_mg)
conditional_effects(fit_mg, effects = "treatment")conditional_effects(fit_mg, effects = "date:treatment")emm_treatments <- emmeans(fit_mg, ~ treatment, type = "response")NOTE: Results may be misleading due to involvement in interactions
emm_treatments treatment response lower.HPD upper.HPD
MC0 51.7 31.1 77.2
MC1 77.3 46.0 114.8
MC2 97.6 59.8 145.5
PC0 32.2 19.8 48.3
PC1 35.7 21.5 52.9
PC2 35.4 21.2 52.9
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
pairs(emm_treatments) contrast ratio lower.HPD upper.HPD
MC0 / MC1 0.670 0.433 0.967
MC0 / MC2 0.530 0.336 0.755
MC0 / PC0 1.608 1.039 2.333
MC0 / PC1 1.447 0.930 2.082
MC0 / PC2 1.462 0.932 2.091
MC1 / MC2 0.791 0.512 1.130
MC1 / PC0 2.407 1.534 3.463
MC1 / PC1 2.161 1.391 3.082
MC1 / PC2 2.176 1.375 3.086
MC2 / PC0 3.036 1.931 4.361
MC2 / PC1 2.728 1.766 3.907
MC2 / PC2 2.757 1.785 3.921
PC0 / PC1 0.901 0.581 1.297
PC0 / PC2 0.906 0.578 1.294
PC1 / PC2 1.006 0.652 1.443
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
plot(emm_treatments) + theme_bw()emmip(fit_mg, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data,
aes(x = date, y = mg, color = as.factor(mg_censored)),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
theme_bw()Warning: Removed 87 rows containing missing values or values outside the scale range
(`geom_point()`).
fit_glmer_gamma_mg <- glmer(
mg ~ date * treatment +
(1 | site_block_id) +
(1 | trt_site_block_id),
data = lc_model_data %>% filter(study_year == 1),
family = Gamma(link = "log"),
# control = glmerControl(optimizer = "bobyqa"), # default optimizer did not converge
na.action = na.exclude
)Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
failure to converge in 10000 evaluations
Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0499573 (tol = 0.002, component 1)
summary(fit_glmer_gamma_mg)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Gamma ( log )
Formula: mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data %>% filter(study_year == 1)
AIC BIC logLik deviance df.resid
3664.3 3891.8 -1775.1 3550.3 343
Scaled residuals:
Min 1Q Median 3Q Max
-2.8607 -0.5774 -0.0216 0.4681 3.8763
Random effects:
Groups Name Variance Std.Dev.
trt_site_block_id (Intercept) 0.04994 0.2235
site_block_id (Intercept) 0.02173 0.1474
Residual 0.12052 0.3472
Number of obs: 400, groups: trt_site_block_id, 48; site_block_id, 8
Fixed effects:
Estimate Std. Error t value Pr(>|z|)
(Intercept) 3.91105 0.25136 15.560 < 2e-16 ***
date2024-01-12 0.21160 0.20181 1.049 0.294392
date2024-01-23 0.22794 0.19553 1.166 0.243703
date2024-02-07 -0.03191 0.19650 -0.162 0.870998
date2024-02-21 0.07750 0.20333 0.381 0.703102
date2024-03-06 0.19660 0.20321 0.967 0.333314
date2024-03-22 0.39571 0.20300 1.949 0.051257 .
date2024-04-04 0.46917 0.21143 2.219 0.026480 *
date2024-04-16 0.28751 0.20094 1.431 0.152480
treatmentMC1 0.81928 0.29224 2.803 0.005056 **
treatmentMC2 0.54335 0.29847 1.820 0.068692 .
treatmentPC0 -0.08276 0.29335 -0.282 0.777855
treatmentPC1 0.37059 0.29228 1.268 0.204831
treatmentPC2 -0.14076 0.31568 -0.446 0.655674
date2024-01-12:treatmentMC1 -0.27861 0.27063 -1.029 0.303254
date2024-01-23:treatmentMC1 -0.27466 0.26672 -1.030 0.303118
date2024-02-07:treatmentMC1 -0.54121 0.26338 -2.055 0.039889 *
date2024-02-21:treatmentMC1 -0.69731 0.26914 -2.591 0.009572 **
date2024-03-06:treatmentMC1 -0.81882 0.26931 -3.040 0.002362 **
date2024-03-22:treatmentMC1 -0.82064 0.26925 -3.048 0.002305 **
date2024-04-04:treatmentMC1 -0.87750 0.28620 -3.066 0.002169 **
date2024-04-16:treatmentMC1 -0.97674 0.26785 -3.647 0.000266 ***
date2024-01-12:treatmentMC2 -0.11081 0.27352 -0.405 0.685374
date2024-01-23:treatmentMC2 0.09168 0.26912 0.341 0.733342
date2024-02-07:treatmentMC2 -0.17715 0.27082 -0.654 0.513026
date2024-02-21:treatmentMC2 -0.10351 0.27572 -0.375 0.707338
date2024-03-06:treatmentMC2 -0.10454 0.27602 -0.379 0.704894
date2024-03-22:treatmentMC2 -0.12819 0.27576 -0.465 0.642026
date2024-04-04:treatmentMC2 -0.17171 0.28198 -0.609 0.542548
date2024-04-16:treatmentMC2 -0.29992 0.27906 -1.075 0.282495
date2024-01-12:treatmentPC0 -0.26652 0.26820 -0.994 0.320350
date2024-01-23:treatmentPC0 -0.49111 0.26410 -1.860 0.062950 .
date2024-02-07:treatmentPC0 -0.34041 0.26482 -1.285 0.198634
date2024-02-21:treatmentPC0 -0.28518 0.26958 -1.058 0.290120
date2024-03-06:treatmentPC0 -0.32802 0.26941 -1.218 0.223397
date2024-03-22:treatmentPC0 -0.50855 0.26957 -1.887 0.059222 .
date2024-04-04:treatmentPC0 -0.51908 0.27641 -1.878 0.060388 .
date2024-04-16:treatmentPC0 -0.40947 0.26879 -1.523 0.127658
date2024-01-12:treatmentPC1 -0.44969 0.27120 -1.658 0.097280 .
date2024-01-23:treatmentPC1 -0.54340 0.26740 -2.032 0.042137 *
date2024-02-07:treatmentPC1 -0.83036 0.26432 -3.142 0.001681 **
date2024-02-21:treatmentPC1 -0.71730 0.26926 -2.664 0.007724 **
date2024-03-06:treatmentPC1 -0.65606 0.26952 -2.434 0.014924 *
date2024-03-22:treatmentPC1 -0.73684 0.26996 -2.729 0.006345 **
date2024-04-04:treatmentPC1 -0.81262 0.28175 -2.884 0.003925 **
date2024-04-16:treatmentPC1 -0.91260 0.27643 -3.301 0.000962 ***
date2024-01-12:treatmentPC2 -0.08284 0.29625 -0.280 0.779770
date2024-01-23:treatmentPC2 -0.34560 0.28816 -1.199 0.230409
date2024-02-07:treatmentPC2 -0.52254 0.29101 -1.796 0.072560 .
date2024-02-21:treatmentPC2 -0.35759 0.29502 -1.212 0.225484
date2024-03-06:treatmentPC2 -0.37510 0.29881 -1.255 0.209362
date2024-03-22:treatmentPC2 -0.36350 0.29492 -1.233 0.217737
date2024-04-04:treatmentPC2 -0.38965 0.30512 -1.277 0.201590
date2024-04-16:treatmentPC2 -0.44361 0.29343 -1.512 0.130582
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation matrix not shown by default, as p = 54 > 12.
Use print(value, correlation=TRUE) or
vcov(value) if you need it
optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations)
Model failed to converge with max|grad| = 0.0499573 (tol = 0.002, component 1)
failure to converge in 10000 evaluations
S(fit_mg)Warning: There were 149 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gamma
Links: mu = log; shape = identity
Formula: mg | cens(mg_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data (Number of observations: 633)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~site_block_id (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.40 0.18 0.17 0.87 1.00 2032 743
~trt_site_block_id (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.38 0.05 0.29 0.50 1.00 6068 9749
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 3.90 0.30 3.29 4.45 1.00 869
date2024M01M12 0.20 0.24 -0.28 0.68 1.00 2141
date2024M01M23 0.25 0.24 -0.23 0.71 1.00 1933
date2024M02M07 -0.01 0.24 -0.48 0.46 1.00 2122
date2024M02M21 0.12 0.25 -0.36 0.62 1.00 2153
date2024M03M06 0.19 0.25 -0.30 0.66 1.00 1806
date2024M03M22 0.38 0.24 -0.09 0.86 1.00 1848
date2024M04M04 0.48 0.26 -0.02 0.98 1.00 2203
date2024M04M16 0.34 0.24 -0.14 0.81 1.00 2128
date2025M02M08 -0.66 0.24 -1.14 -0.19 1.00 1980
date2025M02M15 -0.48 0.25 -0.97 0.00 1.00 2085
date2025M03M04 -0.35 0.26 -0.85 0.16 1.00 2299
date2025M03M16 -0.14 0.24 -0.61 0.33 1.00 2104
date2025M03M27 0.03 0.26 -0.47 0.53 1.00 2299
date2025M04M11 0.21 0.27 -0.31 0.75 1.00 1922
treatmentMC1 0.86 0.31 0.25 1.47 1.00 2513
treatmentMC2 0.58 0.32 -0.04 1.22 1.00 2455
treatmentPC0 0.01 0.32 -0.62 0.62 1.00 2577
treatmentPC1 0.48 0.31 -0.12 1.09 1.00 2187
treatmentPC2 0.10 0.35 -0.57 0.79 1.00 1728
date2024M01M12:treatmentMC1 -0.27 0.33 -0.92 0.38 1.00 2988
date2024M01M23:treatmentMC1 -0.32 0.33 -0.96 0.33 1.00 2830
date2024M02M07:treatmentMC1 -0.55 0.32 -1.20 0.07 1.00 2961
date2024M02M21:treatmentMC1 -0.73 0.33 -1.39 -0.09 1.00 2731
date2024M03M06:treatmentMC1 -0.81 0.33 -1.47 -0.15 1.00 2874
date2024M03M22:treatmentMC1 -0.80 0.33 -1.44 -0.15 1.00 2825
date2024M04M04:treatmentMC1 -0.92 0.35 -1.60 -0.23 1.00 2874
date2024M04M16:treatmentMC1 -1.01 0.33 -1.66 -0.37 1.00 2227
date2025M02M08:treatmentMC1 0.01 0.34 -0.67 0.66 1.00 2500
date2025M02M15:treatmentMC1 -0.07 0.33 -0.72 0.57 1.00 2770
date2025M03M04:treatmentMC1 -0.10 0.34 -0.78 0.56 1.00 3144
date2025M03M16:treatmentMC1 -0.47 0.32 -1.12 0.16 1.00 2787
date2025M03M27:treatmentMC1 -0.50 0.34 -1.19 0.16 1.00 2750
date2025M04M11:treatmentMC1 -0.40 0.36 -1.11 0.29 1.00 3045
date2024M01M12:treatmentMC2 -0.04 0.33 -0.69 0.62 1.00 2875
date2024M01M23:treatmentMC2 0.06 0.33 -0.59 0.70 1.00 2775
date2024M02M07:treatmentMC2 -0.19 0.33 -0.84 0.45 1.00 2798
date2024M02M21:treatmentMC2 -0.14 0.34 -0.82 0.51 1.00 2710
date2024M03M06:treatmentMC2 -0.10 0.34 -0.76 0.57 1.00 2825
date2024M03M22:treatmentMC2 -0.11 0.33 -0.77 0.53 1.00 2511
date2024M04M04:treatmentMC2 -0.19 0.34 -0.86 0.48 1.00 2886
date2024M04M16:treatmentMC2 -0.37 0.34 -1.04 0.29 1.00 2516
date2025M02M08:treatmentMC2 0.19 0.35 -0.49 0.86 1.00 2463
date2025M02M15:treatmentMC2 0.66 0.34 -0.01 1.31 1.00 2887
date2025M03M04:treatmentMC2 0.55 0.35 -0.15 1.23 1.00 3094
date2025M03M16:treatmentMC2 0.21 0.33 -0.44 0.85 1.00 2700
date2025M03M27:treatmentMC2 0.12 0.34 -0.56 0.79 1.00 3062
date2025M04M11:treatmentMC2 0.12 0.39 -0.64 0.88 1.00 3304
date2024M01M12:treatmentPC0 -0.31 0.33 -0.95 0.35 1.00 2992
date2024M01M23:treatmentPC0 -0.56 0.32 -1.19 0.07 1.00 2745
date2024M02M07:treatmentPC0 -0.40 0.33 -1.05 0.23 1.00 2922
date2024M02M21:treatmentPC0 -0.37 0.33 -1.02 0.28 1.00 3116
date2024M03M06:treatmentPC0 -0.36 0.33 -1.00 0.30 1.00 2483
date2024M03M22:treatmentPC0 -0.54 0.33 -1.19 0.11 1.00 2706
date2024M04M04:treatmentPC0 -0.59 0.34 -1.25 0.09 1.00 2953
date2024M04M16:treatmentPC0 -0.51 0.33 -1.17 0.13 1.00 3257
date2025M02M08:treatmentPC0 -0.59 0.39 -1.35 0.17 1.00 3075
date2025M02M15:treatmentPC0 -0.55 0.33 -1.21 0.10 1.00 2969
date2025M03M04:treatmentPC0 -0.46 0.36 -1.17 0.26 1.00 3014
date2025M03M16:treatmentPC0 -0.81 0.34 -1.47 -0.14 1.00 3140
date2025M03M27:treatmentPC0 -0.67 0.38 -1.39 0.08 1.00 3335
date2025M04M11:treatmentPC0 -0.53 0.37 -1.25 0.19 1.00 2247
date2024M01M12:treatmentPC1 -0.46 0.33 -1.11 0.19 1.00 2995
date2024M01M23:treatmentPC1 -0.60 0.33 -1.23 0.05 1.00 2892
date2024M02M07:treatmentPC1 -0.91 0.32 -1.54 -0.29 1.00 2752
date2024M02M21:treatmentPC1 -0.84 0.33 -1.48 -0.21 1.00 2476
date2024M03M06:treatmentPC1 -0.73 0.33 -1.37 -0.08 1.00 2880
date2024M03M22:treatmentPC1 -0.82 0.33 -1.46 -0.18 1.00 2772
date2024M04M04:treatmentPC1 -0.95 0.34 -1.61 -0.28 1.00 2979
date2024M04M16:treatmentPC1 -1.03 0.34 -1.70 -0.37 1.00 2945
date2025M02M08:treatmentPC1 -1.01 0.33 -1.66 -0.36 1.00 2785
date2025M02M15:treatmentPC1 -1.05 0.34 -1.72 -0.39 1.00 2882
date2025M03M04:treatmentPC1 -1.12 0.34 -1.80 -0.45 1.00 3182
date2025M03M16:treatmentPC1 -1.05 0.34 -1.71 -0.39 1.00 3087
date2025M03M27:treatmentPC1 -1.04 0.35 -1.72 -0.34 1.00 3212
date2025M04M11:treatmentPC1 -1.23 0.40 -2.00 -0.44 1.00 3125
date2024M01M12:treatmentPC2 -0.29 0.36 -1.00 0.43 1.00 2778
date2024M01M23:treatmentPC2 -0.55 0.36 -1.25 0.14 1.00 2333
date2024M02M07:treatmentPC2 -0.78 0.36 -1.47 -0.07 1.00 2159
date2024M02M21:treatmentPC2 -0.61 0.36 -1.31 0.11 1.00 2491
date2024M03M06:treatmentPC2 -0.58 0.37 -1.30 0.13 1.00 1856
date2024M03M22:treatmentPC2 -0.55 0.36 -1.26 0.15 1.00 2027
date2024M04M04:treatmentPC2 -0.62 0.37 -1.34 0.12 1.00 2706
date2024M04M16:treatmentPC2 -0.68 0.36 -1.39 0.02 1.00 3141
date2025M02M08:treatmentPC2 -0.65 0.37 -1.37 0.07 1.00 2916
date2025M02M15:treatmentPC2 -0.36 0.36 -1.07 0.35 1.00 2638
date2025M03M04:treatmentPC2 -0.22 0.38 -0.96 0.52 1.00 2541
date2025M03M16:treatmentPC2 -0.40 0.36 -1.11 0.30 1.00 2830
date2025M03M27:treatmentPC2 -0.38 0.38 -1.14 0.36 1.00 2334
date2025M04M11:treatmentPC2 -0.56 0.43 -1.38 0.30 1.00 2570
Tail_ESS
Intercept 489
date2024M01M12 4611
date2024M01M23 4161
date2024M02M07 4299
date2024M02M21 5113
date2024M03M06 5267
date2024M03M22 5165
date2024M04M04 5257
date2024M04M16 5252
date2025M02M08 4562
date2025M02M15 5271
date2025M03M04 5281
date2025M03M16 5087
date2025M03M27 5755
date2025M04M11 5094
treatmentMC1 5223
treatmentMC2 6142
treatmentPC0 6150
treatmentPC1 5902
treatmentPC2 1876
date2024M01M12:treatmentMC1 7291
date2024M01M23:treatmentMC1 6716
date2024M02M07:treatmentMC1 6611
date2024M02M21:treatmentMC1 5988
date2024M03M06:treatmentMC1 6631
date2024M03M22:treatmentMC1 6946
date2024M04M04:treatmentMC1 6233
date2024M04M16:treatmentMC1 5203
date2025M02M08:treatmentMC1 6029
date2025M02M15:treatmentMC1 5923
date2025M03M04:treatmentMC1 6228
date2025M03M16:treatmentMC1 5853
date2025M03M27:treatmentMC1 7492
date2025M04M11:treatmentMC1 7125
date2024M01M12:treatmentMC2 7081
date2024M01M23:treatmentMC2 6424
date2024M02M07:treatmentMC2 6043
date2024M02M21:treatmentMC2 6387
date2024M03M06:treatmentMC2 7069
date2024M03M22:treatmentMC2 5811
date2024M04M04:treatmentMC2 6399
date2024M04M16:treatmentMC2 7370
date2025M02M08:treatmentMC2 6477
date2025M02M15:treatmentMC2 6186
date2025M03M04:treatmentMC2 6360
date2025M03M16:treatmentMC2 6414
date2025M03M27:treatmentMC2 7132
date2025M04M11:treatmentMC2 6804
date2024M01M12:treatmentPC0 6866
date2024M01M23:treatmentPC0 7129
date2024M02M07:treatmentPC0 6789
date2024M02M21:treatmentPC0 6619
date2024M03M06:treatmentPC0 5943
date2024M03M22:treatmentPC0 7353
date2024M04M04:treatmentPC0 6545
date2024M04M16:treatmentPC0 7638
date2025M02M08:treatmentPC0 8050
date2025M02M15:treatmentPC0 6781
date2025M03M04:treatmentPC0 7508
date2025M03M16:treatmentPC0 6850
date2025M03M27:treatmentPC0 7941
date2025M04M11:treatmentPC0 5367
date2024M01M12:treatmentPC1 6645
date2024M01M23:treatmentPC1 6957
date2024M02M07:treatmentPC1 7158
date2024M02M21:treatmentPC1 6270
date2024M03M06:treatmentPC1 7322
date2024M03M22:treatmentPC1 6564
date2024M04M04:treatmentPC1 6755
date2024M04M16:treatmentPC1 7499
date2025M02M08:treatmentPC1 6374
date2025M02M15:treatmentPC1 7465
date2025M03M04:treatmentPC1 7572
date2025M03M16:treatmentPC1 7202
date2025M03M27:treatmentPC1 7874
date2025M04M11:treatmentPC1 8764
date2024M01M12:treatmentPC2 6428
date2024M01M23:treatmentPC2 6396
date2024M02M07:treatmentPC2 6143
date2024M02M21:treatmentPC2 6478
date2024M03M06:treatmentPC2 3113
date2024M03M22:treatmentPC2 4921
date2024M04M04:treatmentPC2 7305
date2024M04M16:treatmentPC2 6883
date2025M02M08:treatmentPC2 7278
date2025M02M15:treatmentPC2 7244
date2025M03M04:treatmentPC2 6715
date2025M03M16:treatmentPC2 6105
date2025M03M27:treatmentPC2 3233
date2025M04M11:treatmentPC2 7080
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 5.78 0.35 5.11 6.48 1.00 13355 11217
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fitted(fit_glmer_gamma_mg), residuals(fit_glmer_gamma_mg),
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residuals vs. Fitted Values")
abline(h = 0, col = "red", lty = 2)emmip(fit_glmer_gamma_mg, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data %>% filter(study_year == 1),
aes(x = date, y = mg, color = treatment),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
theme_bw()Warning: Removed 32 rows containing missing values or values outside the scale range
(`geom_point()`).
emmip(fit_mg, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data %>% filter(study_year == 1),
aes(x = date, y = mg, color = treatment),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
theme_bw()Warning: Removed 32 rows containing missing values or values outside the scale range
(`geom_point()`).
summary(fit_p)Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
careful when analysing the results! We recommend running more iterations and/or
setting stronger priors.
Warning: There were 3210 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gamma
Links: mu = log; shape = identity
Formula: p | cens(p_censored) ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Data: lc_model_data (Number of observations: 626)
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~site_block_id (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.15 0.12 0.01 0.46 1.01 460 632
~trt_site_block_id (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.59 0.10 0.42 0.83 1.00 534 947
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept -1.23 0.59 -2.34 -0.13 1.09 38
date2024M01M12 -0.51 0.71 -1.89 0.84 1.10 36
date2024M01M23 -1.28 0.71 -2.68 0.03 1.09 46
date2024M02M07 -2.05 0.70 -3.39 -0.74 1.08 44
date2024M02M21 -2.69 0.73 -4.07 -1.29 1.11 32
date2024M03M06 -2.22 0.73 -3.69 -0.83 1.06 53
date2024M03M22 -1.31 0.73 -2.73 0.20 1.08 41
date2024M04M04 -1.37 0.74 -2.77 0.11 1.08 50
date2024M04M16 -2.19 0.71 -3.51 -0.79 1.07 51
date2025M02M08 -0.13 0.70 -1.49 1.20 1.07 55
date2025M02M15 -1.31 0.75 -2.80 0.13 1.06 59
date2025M03M04 -284.60 205.13 -678.11 -11.41 1.02 117
date2025M03M16 -1.42 0.76 -2.95 0.06 1.09 36
date2025M03M27 -0.92 0.77 -2.36 0.61 1.06 58
date2025M04M11 -0.07 0.87 -1.74 1.70 1.07 49
treatmentMC1 -0.15 0.75 -1.59 1.31 1.09 37
treatmentMC2 -0.36 0.74 -1.80 1.13 1.09 37
treatmentPC0 0.54 0.80 -0.93 2.21 1.06 54
treatmentPC1 0.14 0.76 -1.32 1.61 1.06 57
treatmentPC2 0.18 0.93 -1.43 2.20 1.05 62
date2024M01M12:treatmentMC1 -0.25 0.95 -2.12 1.65 1.09 39
date2024M01M23:treatmentMC1 0.07 0.95 -1.78 1.92 1.08 51
date2024M02M07:treatmentMC1 -0.36 0.94 -2.20 1.54 1.07 48
date2024M02M21:treatmentMC1 0.37 0.96 -1.53 2.22 1.10 33
date2024M03M06:treatmentMC1 -0.49 0.95 -2.31 1.46 1.07 54
date2024M03M22:treatmentMC1 -1.74 0.98 -3.65 0.24 1.08 36
date2024M04M04:treatmentMC1 -1.41 0.96 -3.30 0.48 1.08 50
date2024M04M16:treatmentMC1 -0.62 0.94 -2.53 1.14 1.06 53
date2025M02M08:treatmentMC1 0.32 0.98 -1.58 2.16 1.06 68
date2025M02M15:treatmentMC1 -0.38 1.00 -2.26 1.58 1.05 65
date2025M03M04:treatmentMC1 282.26 205.12 8.99 675.95 1.02 117
date2025M03M16:treatmentMC1 -233.74 186.65 -650.85 -8.95 1.00 490
date2025M03M27:treatmentMC1 -235.65 186.00 -654.17 -9.65 1.00 496
date2025M04M11:treatmentMC1 -0.44 1.14 -2.65 1.75 1.06 55
date2024M01M12:treatmentMC2 -0.03 0.93 -1.81 1.80 1.07 56
date2024M01M23:treatmentMC2 -0.37 0.92 -2.15 1.41 1.07 59
date2024M02M07:treatmentMC2 -0.57 0.92 -2.40 1.21 1.07 58
date2024M02M21:treatmentMC2 0.19 0.91 -1.63 1.97 1.08 43
date2024M03M06:treatmentMC2 -0.53 0.96 -2.39 1.32 1.05 72
date2024M03M22:treatmentMC2 -1.96 0.93 -3.87 -0.18 1.07 49
date2024M04M04:treatmentMC2 -1.18 0.93 -3.01 0.63 1.07 58
date2024M04M16:treatmentMC2 0.64 0.96 -1.20 2.52 1.06 60
date2025M02M08:treatmentMC2 -0.87 0.94 -2.70 0.94 1.05 75
date2025M02M15:treatmentMC2 0.14 0.97 -1.73 2.04 1.05 67
date2025M03M04:treatmentMC2 31.68 275.55 -512.93 561.27 1.01 201
date2025M03M16:treatmentMC2 -230.00 181.27 -644.10 -8.83 1.00 572
date2025M03M27:treatmentMC2 -0.13 1.01 -2.16 1.79 1.05 68
date2025M04M11:treatmentMC2 -0.96 1.20 -3.26 1.46 1.04 85
date2024M01M12:treatmentPC0 -1.08 0.97 -3.03 0.77 1.05 69
date2024M01M23:treatmentPC0 -1.33 0.97 -3.25 0.50 1.05 65
date2024M02M07:treatmentPC0 -1.33 0.99 -3.21 0.56 1.04 76
date2024M02M21:treatmentPC0 -0.00 1.01 -1.98 1.87 1.05 75
date2024M03M06:treatmentPC0 -0.97 1.02 -2.98 1.04 1.04 67
date2024M03M22:treatmentPC0 -2.23 1.00 -4.25 -0.36 1.05 74
date2024M04M04:treatmentPC0 -1.29 1.00 -3.26 0.59 1.05 77
date2024M04M16:treatmentPC0 -0.25 0.97 -2.15 1.58 1.04 76
date2025M02M08:treatmentPC0 -0.11 1.20 -2.42 2.23 1.03 96
date2025M02M15:treatmentPC0 0.61 1.04 -1.44 2.56 1.04 67
date2025M03M04:treatmentPC0 283.09 205.16 9.75 676.78 1.02 117
date2025M03M16:treatmentPC0 -258.26 191.44 -658.91 -10.67 1.00 775
date2025M03M27:treatmentPC0 -0.10 1.15 -2.34 2.17 1.03 115
date2025M04M11:treatmentPC0 -0.06 1.22 -2.49 2.27 1.04 81
date2024M01M12:treatmentPC1 -0.27 0.93 -2.05 1.54 1.06 70
date2024M01M23:treatmentPC1 -0.51 0.94 -2.33 1.39 1.06 66
date2024M02M07:treatmentPC1 0.01 0.90 -1.75 1.72 1.05 73
date2024M02M21:treatmentPC1 0.27 0.93 -1.54 2.10 1.07 60
date2024M03M06:treatmentPC1 -0.37 0.94 -2.14 1.47 1.05 68
date2024M03M22:treatmentPC1 -1.56 0.94 -3.41 0.28 1.05 82
date2024M04M04:treatmentPC1 -0.95 0.96 -2.82 0.94 1.05 86
date2024M04M16:treatmentPC1 1.19 0.97 -0.73 3.07 1.04 90
date2025M02M08:treatmentPC1 1.25 0.97 -0.63 3.17 1.04 89
date2025M02M15:treatmentPC1 1.00 0.98 -0.93 2.93 1.03 86
date2025M03M04:treatmentPC1 283.12 205.13 9.83 676.79 1.02 117
date2025M03M16:treatmentPC1 0.64 0.99 -1.30 2.60 1.06 79
date2025M03M27:treatmentPC1 0.74 1.04 -1.32 2.75 1.04 92
date2025M04M11:treatmentPC1 1.19 1.23 -1.19 3.61 1.03 99
date2024M01M12:treatmentPC2 -0.25 1.07 -2.43 1.71 1.06 55
date2024M01M23:treatmentPC2 0.11 1.06 -2.12 2.04 1.06 53
date2024M02M07:treatmentPC2 0.33 1.06 -1.88 2.22 1.06 58
date2024M02M21:treatmentPC2 -0.03 1.06 -2.13 1.93 1.06 47
date2024M03M06:treatmentPC2 -0.33 1.06 -2.53 1.64 1.04 69
date2024M03M22:treatmentPC2 -0.63 1.08 -2.88 1.35 1.05 58
date2024M04M04:treatmentPC2 0.25 1.07 -1.98 2.27 1.05 71
date2024M04M16:treatmentPC2 0.65 1.08 -1.56 2.63 1.05 71
date2025M02M08:treatmentPC2 0.71 1.08 -1.47 2.79 1.05 69
date2025M02M15:treatmentPC2 1.23 1.07 -0.91 3.24 1.04 78
date2025M03M04:treatmentPC2 283.98 205.09 10.88 677.53 1.02 117
date2025M03M16:treatmentPC2 1.42 1.09 -0.76 3.51 1.05 56
date2025M03M27:treatmentPC2 1.38 1.16 -1.04 3.58 1.04 74
date2025M04M11:treatmentPC2 1.30 1.33 -1.33 3.90 1.04 74
Tail_ESS
Intercept 152
date2024M01M12 166
date2024M01M23 253
date2024M02M07 192
date2024M02M21 351
date2024M03M06 96
date2024M03M22 327
date2024M04M04 243
date2024M04M16 184
date2025M02M08 312
date2025M02M15 106
date2025M03M04 366
date2025M03M16 100
date2025M03M27 558
date2025M04M11 135
treatmentMC1 231
treatmentMC2 184
treatmentPC0 194
treatmentPC1 299
treatmentPC2 154
date2024M01M12:treatmentMC1 202
date2024M01M23:treatmentMC1 282
date2024M02M07:treatmentMC1 223
date2024M02M21:treatmentMC1 317
date2024M03M06:treatmentMC1 129
date2024M03M22:treatmentMC1 338
date2024M04M04:treatmentMC1 264
date2024M04M16:treatmentMC1 240
date2025M02M08:treatmentMC1 385
date2025M02M15:treatmentMC1 211
date2025M03M04:treatmentMC1 361
date2025M03M16:treatmentMC1 1068
date2025M03M27:treatmentMC1 956
date2025M04M11:treatmentMC1 157
date2024M01M12:treatmentMC2 346
date2024M01M23:treatmentMC2 303
date2024M02M07:treatmentMC2 248
date2024M02M21:treatmentMC2 216
date2024M03M06:treatmentMC2 176
date2024M03M22:treatmentMC2 316
date2024M04M04:treatmentMC2 511
date2024M04M16:treatmentMC2 595
date2025M02M08:treatmentMC2 377
date2025M02M15:treatmentMC2 210
date2025M03M04:treatmentMC2 670
date2025M03M16:treatmentMC2 993
date2025M03M27:treatmentMC2 590
date2025M04M11:treatmentMC2 579
date2024M01M12:treatmentPC0 310
date2024M01M23:treatmentPC0 232
date2024M02M07:treatmentPC0 280
date2024M02M21:treatmentPC0 407
date2024M03M06:treatmentPC0 186
date2024M03M22:treatmentPC0 281
date2024M04M04:treatmentPC0 314
date2024M04M16:treatmentPC0 183
date2025M02M08:treatmentPC0 471
date2025M02M15:treatmentPC0 254
date2025M03M04:treatmentPC0 366
date2025M03M16:treatmentPC0 1263
date2025M03M27:treatmentPC0 654
date2025M04M11:treatmentPC0 478
date2024M01M12:treatmentPC1 413
date2024M01M23:treatmentPC1 323
date2024M02M07:treatmentPC1 462
date2024M02M21:treatmentPC1 374
date2024M03M06:treatmentPC1 213
date2024M03M22:treatmentPC1 400
date2024M04M04:treatmentPC1 388
date2024M04M16:treatmentPC1 405
date2025M02M08:treatmentPC1 365
date2025M02M15:treatmentPC1 265
date2025M03M04:treatmentPC1 348
date2025M03M16:treatmentPC1 560
date2025M03M27:treatmentPC1 679
date2025M04M11:treatmentPC1 559
date2024M01M12:treatmentPC2 177
date2024M01M23:treatmentPC2 203
date2024M02M07:treatmentPC2 168
date2024M02M21:treatmentPC2 186
date2024M03M06:treatmentPC2 231
date2024M03M22:treatmentPC2 226
date2024M04M04:treatmentPC2 315
date2024M04M16:treatmentPC2 260
date2025M02M08:treatmentPC2 290
date2025M02M15:treatmentPC2 205
date2025M03M04:treatmentPC2 363
date2025M03M16:treatmentPC2 246
date2025M03M27:treatmentPC2 212
date2025M04M11:treatmentPC2 612
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.80 0.06 0.69 0.91 1.00 863 1421
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
# neff_ratio(fit_p)
# prior_summary(fit_p)
# plot(fit, type = "trace", nvariables = 1)
pp_check(fit_p) + xlim(0, 1)Using 10 posterior draws for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 189 rows containing non-finite outside the scale range
(`stat_density()`).
Warning: Removed 18 rows containing non-finite outside the scale range
(`stat_density()`).
pp_check(fit_p, type = "dens_overlay", ndraws = 100) + xlim(0, 1)Warning: Censored responses are not shown in 'pp_check'.
Warning: Removed 1818 rows containing non-finite outside the scale range
(`stat_density()`).
Warning: Removed 18 rows containing non-finite outside the scale range
(`stat_density()`).
# loo(fit_p)
conditional_effects(fit_p, effects = "treatment")conditional_effects(fit_p, effects = "date:treatment")emm_treatments <- emmeans(fit_p, ~ treatment, type="response")NOTE: Results may be misleading due to involvement in interactions
emm_treatments treatment response lower.HPD upper.HPD
MC0 1.00e-08 0.0000 2.37e-02
MC1 0.00e+00 0.0000 8.72e-05
MC2 0.00e+00 0.0000 4.81e-05
PC0 5.00e-08 0.0000 2.28e-02
PC1 1.12e-01 0.0620 1.76e-01
PC2 1.57e-01 0.0831 2.54e-01
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
pairs(emm_treatments) contrast ratio lower.HPD upper.HPD
MC0 / MC1 146514.1 0.000 2.59e+22
MC0 / MC2 463834.4 0.000 1.79e+22
MC0 / PC0 0.3 0.000 3.82e+12
MC0 / PC1 0.0 0.000 0.00e+00
MC0 / PC2 0.0 0.000 0.00e+00
MC1 / MC2 3.0 0.000 4.38e+17
MC1 / PC0 0.0 0.000 2.34e+08
MC1 / PC1 0.0 0.000 0.00e+00
MC1 / PC2 0.0 0.000 0.00e+00
MC2 / PC0 0.0 0.000 7.10e+07
MC2 / PC1 0.0 0.000 0.00e+00
MC2 / PC2 0.0 0.000 0.00e+00
PC0 / PC1 0.0 0.000 0.00e+00
PC0 / PC2 0.0 0.000 0.00e+00
PC1 / PC2 0.7 0.306 1.00e+00
Results are averaged over the levels of: date
Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95
plot(emm_treatments) + theme_bw()emmip(fit_p, ~ date | treatment, type = "response", CIs = T) +
geom_point(data = lc_model_data,
aes(x = date, y = p, color = as.factor(p_censored)),
alpha = 0.7,
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.1)) +
# ylim(0, 4) +
theme_bw()Warning: Removed 94 rows containing missing values or values outside the scale range
(`geom_point()`).
fit_no3 <- glmmTMB(
no3 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
zi = ~ treatment,
family = glmmTMB::lognormal(link = "log"),
data = df_model
)Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
S(fit_no3) Family: lognormal ( log )
Formula:
no3 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation: ~treatment
Data: df_model
AIC BIC logLik -2*log(L) df.resid
NA NA NA NA 343
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site_block_id (Intercept) 0.18431 0.4293
trt_site_block_id (Intercept) 0.07486 0.2736
Number of obs: 406, groups: site_block_id, 8; trt_site_block_id, 48
Dispersion parameter for lognormal family (): 43.2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.94296 NaN NaN NaN
date2024-01-12 -0.46589 NaN NaN NaN
date2024-01-23 -2.24045 NaN NaN NaN
date2024-02-07 -1.16452 NaN NaN NaN
date2024-02-21 -0.93779 NaN NaN NaN
date2024-03-06 -0.58716 NaN NaN NaN
date2024-03-22 -2.06985 NaN NaN NaN
date2024-04-04 -0.85750 NaN NaN NaN
date2024-04-16 -0.35219 NaN NaN NaN
treatmentMC1 -0.95229 NaN NaN NaN
treatmentMC2 -0.18415 NaN NaN NaN
treatmentPC0 -0.59872 NaN NaN NaN
treatmentPC1 -0.43191 NaN NaN NaN
treatmentPC2 -0.26465 NaN NaN NaN
date2024-01-12:treatmentMC1 0.35628 NaN NaN NaN
date2024-01-23:treatmentMC1 2.13208 NaN NaN NaN
date2024-02-07:treatmentMC1 -0.03654 NaN NaN NaN
date2024-02-21:treatmentMC1 0.00000 NaN NaN NaN
date2024-03-06:treatmentMC1 0.00000 NaN NaN NaN
date2024-03-22:treatmentMC1 0.71224 NaN NaN NaN
date2024-04-04:treatmentMC1 0.00000 NaN NaN NaN
date2024-04-16:treatmentMC1 0.00000 NaN NaN NaN
date2024-01-12:treatmentMC2 -0.51072 NaN NaN NaN
date2024-01-23:treatmentMC2 0.20048 NaN NaN NaN
date2024-02-07:treatmentMC2 0.48826 NaN NaN NaN
date2024-02-21:treatmentMC2 -0.29661 NaN NaN NaN
date2024-03-06:treatmentMC2 -1.29413 NaN NaN NaN
date2024-03-22:treatmentMC2 0.00000 NaN NaN NaN
date2024-04-04:treatmentMC2 0.00000 NaN NaN NaN
date2024-04-16:treatmentMC2 -1.03794 NaN NaN NaN
date2024-01-12:treatmentPC0 -0.12570 NaN NaN NaN
date2024-01-23:treatmentPC0 0.84651 NaN NaN NaN
date2024-02-07:treatmentPC0 0.00000 NaN NaN NaN
date2024-02-21:treatmentPC0 -1.56917 NaN NaN NaN
date2024-03-06:treatmentPC0 0.00000 NaN NaN NaN
date2024-03-22:treatmentPC0 0.00000 NaN NaN NaN
date2024-04-04:treatmentPC0 0.00000 NaN NaN NaN
date2024-04-16:treatmentPC0 0.00000 NaN NaN NaN
date2024-01-12:treatmentPC1 -0.32762 NaN NaN NaN
date2024-01-23:treatmentPC1 0.64563 NaN NaN NaN
date2024-02-07:treatmentPC1 1.88958 NaN NaN NaN
date2024-02-21:treatmentPC1 0.92799 NaN NaN NaN
date2024-03-06:treatmentPC1 1.48469 NaN NaN NaN
date2024-03-22:treatmentPC1 1.94375 NaN NaN NaN
date2024-04-04:treatmentPC1 0.00000 NaN NaN NaN
date2024-04-16:treatmentPC1 0.68575 NaN NaN NaN
date2024-01-12:treatmentPC2 -0.11587 NaN NaN NaN
date2024-01-23:treatmentPC2 0.56902 NaN NaN NaN
date2024-02-07:treatmentPC2 0.00000 NaN NaN NaN
date2024-02-21:treatmentPC2 0.00000 NaN NaN NaN
date2024-03-06:treatmentPC2 -0.77773 NaN NaN NaN
date2024-03-22:treatmentPC2 -0.03610 NaN NaN NaN
date2024-04-04:treatmentPC2 -0.85750 NaN NaN NaN
date2024-04-16:treatmentPC2 0.00000 NaN NaN NaN
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6702 NaN NaN NaN
treatmentMC1 0.1565 NaN NaN NaN
treatmentMC2 0.2974 NaN NaN NaN
treatmentPC0 0.3366 NaN NaN NaN
treatmentPC1 0.3313 NaN NaN NaN
treatmentPC2 0.2141 NaN NaN NaN
simulationOutput <- simulateResiduals(fittedModel = fit_no3, n = 500)
plot(simulationOutput)testUniformity(simulationOutput)
Asymptotic one-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.035847, p-value = 0.6739
alternative hypothesis: two-sided
testOutliers(simulationOutput)
DHARMa outlier test based on exact binomial test with approximate
expectations
data: simulationOutput
outliers at both margin(s) = 2, observations = 406, p-value = 0.6793
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
0.0005971323 0.0176806666
sample estimates:
frequency of outliers (expected: 0.00399201596806387 )
0.004926108
testDispersion(simulationOutput)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.84685, p-value = 1
alternative hypothesis: two.sided
testZeroInflation(simulationOutput)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 0.99872, p-value = 1
alternative hypothesis: two.sided
library(ggeffects)Warning: package 'ggeffects' was built under R version 4.4.3
library(RColorBrewer)
pred_effects <- ggpredict(fit_no3, terms = c("date", "treatment"))You are calculating adjusted predictions on the population-level (i.e.
`type = "fixed"`) for a *generalized* linear mixed model.
This may produce biased estimates due to Jensen's inequality. Consider
setting `bias_correction = TRUE` to correct for this bias.
See also the documentation of the `bias_correction` argument.
plot(pred_effects, colors = "flat") +
theme_bw() +
labs(
x = "Date",
y = "Nitrate",
title = "Predicted Values of Nitrate",
color = "Treatment"
)print(pred_effects)# Predicted values of no3
treatment: MC0
date | Predicted
----------------------
2023-12-27 | 140.18
2024-01-12 | 87.98
2024-01-23 | 14.92
2024-02-07 | 43.75
2024-02-21 | 54.88
2024-03-06 | 77.93
2024-03-22 | 17.69
2024-04-04 | 59.47
2024-04-16 | 98.57
treatment: MC1
date | Predicted
----------------------
2023-12-27 | 54.09
2024-01-12 | 48.48
2024-01-23 | 48.54
2024-02-07 | 16.27
2024-02-21 | 21.18
2024-03-06 | 30.07
2024-03-22 | 13.92
2024-04-04 | 22.95
2024-04-16 | 38.03
treatment: MC2
date | Predicted
----------------------
2023-12-27 | 116.61
2024-01-12 | 43.91
2024-01-23 | 15.16
2024-02-07 | 59.30
2024-02-21 | 33.93
2024-03-06 | 17.77
2024-03-22 | 14.72
2024-04-04 | 49.47
2024-04-16 | 29.04
treatment: PC0
date | Predicted
----------------------
2023-12-27 | 77.03
2024-01-12 | 42.63
2024-01-23 | 19.11
2024-02-07 | 24.04
2024-02-21 | 6.28
2024-03-06 | 42.82
2024-03-22 | 9.72
2024-04-04 | 32.68
2024-04-16 | 54.17
treatment: PC1
date | Predicted
----------------------
2023-12-27 | 91.02
2024-01-12 | 41.16
2024-01-23 | 18.47
2024-02-07 | 187.94
2024-02-21 | 90.13
2024-03-06 | 223.32
2024-03-22 | 80.23
2024-04-04 | 38.61
2024-04-16 | 127.05
treatment: PC2
date | Predicted
----------------------
2023-12-27 | 107.59
2024-01-12 | 60.13
2024-01-23 | 20.22
2024-02-07 | 33.58
2024-02-21 | 42.12
2024-03-06 | 27.48
2024-03-22 | 13.10
2024-04-04 | 19.36
2024-04-16 | 75.65
Adjusted for:
* site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)
fit_nh4 <- glmmTMB(
nh4 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
zi = ~ treatment,
family = glmmTMB::lognormal(link = "log"),
data = df_model
)Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
S(fit_nh4) Family: lognormal ( log )
Formula:
nh4 ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation: ~treatment
Data: df_model
AIC BIC logLik -2*log(L) df.resid
NA NA NA NA 349
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site_block_id (Intercept) 4.380e-11 6.618e-06
trt_site_block_id (Intercept) 3.594e-10 1.896e-05
Number of obs: 412, groups: site_block_id, 8; trt_site_block_id, 48
Dispersion parameter for lognormal family (): 1.13
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.543717 NaN NaN NaN
date2024-01-12 -0.606122 NaN NaN NaN
date2024-01-23 -0.068191 NaN NaN NaN
date2024-02-07 -0.119794 NaN NaN NaN
date2024-02-21 -0.365954 NaN NaN NaN
date2024-03-06 -0.386029 NaN NaN NaN
date2024-03-22 -0.702810 NaN NaN NaN
date2024-04-04 -0.346473 NaN NaN NaN
date2024-04-16 -0.444480 NaN NaN NaN
treatmentMC1 -0.186478 NaN NaN NaN
treatmentMC2 -0.300228 NaN NaN NaN
treatmentPC0 -0.438777 NaN NaN NaN
treatmentPC1 -0.017852 NaN NaN NaN
treatmentPC2 -0.568562 NaN NaN NaN
date2024-01-12:treatmentMC1 -0.009409 NaN NaN NaN
date2024-01-23:treatmentMC1 -0.475401 NaN NaN NaN
date2024-02-07:treatmentMC1 0.000000 NaN NaN NaN
date2024-02-21:treatmentMC1 0.000000 NaN NaN NaN
date2024-03-06:treatmentMC1 0.000000 NaN NaN NaN
date2024-03-22:treatmentMC1 0.000000 NaN NaN NaN
date2024-04-04:treatmentMC1 0.331027 NaN NaN NaN
date2024-04-16:treatmentMC1 0.000000 NaN NaN NaN
date2024-01-12:treatmentMC2 0.362826 NaN NaN NaN
date2024-01-23:treatmentMC2 -0.689456 NaN NaN NaN
date2024-02-07:treatmentMC2 0.000000 NaN NaN NaN
date2024-02-21:treatmentMC2 0.000000 NaN NaN NaN
date2024-03-06:treatmentMC2 0.000000 NaN NaN NaN
date2024-03-22:treatmentMC2 0.000000 NaN NaN NaN
date2024-04-04:treatmentMC2 0.027561 NaN NaN NaN
date2024-04-16:treatmentMC2 0.000000 NaN NaN NaN
date2024-01-12:treatmentPC0 0.000000 NaN NaN NaN
date2024-01-23:treatmentPC0 0.000000 NaN NaN NaN
date2024-02-07:treatmentPC0 0.000000 NaN NaN NaN
date2024-02-21:treatmentPC0 0.000000 NaN NaN NaN
date2024-03-06:treatmentPC0 0.000000 NaN NaN NaN
date2024-03-22:treatmentPC0 0.000000 NaN NaN NaN
date2024-04-04:treatmentPC0 0.026161 NaN NaN NaN
date2024-04-16:treatmentPC0 0.136907 NaN NaN NaN
date2024-01-12:treatmentPC1 -0.251759 NaN NaN NaN
date2024-01-23:treatmentPC1 0.889392 NaN NaN NaN
date2024-02-07:treatmentPC1 0.000000 NaN NaN NaN
date2024-02-21:treatmentPC1 -0.365954 NaN NaN NaN
date2024-03-06:treatmentPC1 0.000000 NaN NaN NaN
date2024-03-22:treatmentPC1 0.000000 NaN NaN NaN
date2024-04-04:treatmentPC1 0.112345 NaN NaN NaN
date2024-04-16:treatmentPC1 0.547386 NaN NaN NaN
date2024-01-12:treatmentPC2 0.275519 NaN NaN NaN
date2024-01-23:treatmentPC2 0.207274 NaN NaN NaN
date2024-02-07:treatmentPC2 -0.119794 NaN NaN NaN
date2024-02-21:treatmentPC2 0.000000 NaN NaN NaN
date2024-03-06:treatmentPC2 0.000000 NaN NaN NaN
date2024-03-22:treatmentPC2 0.000000 NaN NaN NaN
date2024-04-04:treatmentPC2 0.172446 NaN NaN NaN
date2024-04-16:treatmentPC2 0.000000 NaN NaN NaN
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2040 NaN NaN NaN
treatmentMC1 0.2741 NaN NaN NaN
treatmentMC2 0.5710 NaN NaN NaN
treatmentPC0 1.0245 NaN NaN NaN
treatmentPC1 0.2564 NaN NaN NaN
treatmentPC2 0.9445 NaN NaN NaN
simulationOutput <- simulateResiduals(fittedModel = fit_nh4, n = 500)
plot(simulationOutput)testUniformity(simulationOutput)
Asymptotic one-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.033539, p-value = 0.743
alternative hypothesis: two-sided
testOutliers(simulationOutput)
DHARMa outlier test based on exact binomial test with approximate
expectations
data: simulationOutput
outliers at both margin(s) = 2, observations = 412, p-value = 0.6822
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
0.0005884282 0.0174248139
sample estimates:
frequency of outliers (expected: 0.00399201596806387 )
0.004854369
testDispersion(simulationOutput)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.7099, p-value = 0.132
alternative hypothesis: two.sided
testZeroInflation(simulationOutput)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 0.99968, p-value = 1
alternative hypothesis: two.sided
library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_nh4, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") +
theme_bw() +
labs(
x = "Date",
y = "NH4",
title = "Predicted Values of NH4",
color = "Treatment"
)print(pred_effects)# Predicted values of nh4
treatment: MC0
date | Predicted
----------------------
2023-12-27 | 1.72
2024-01-12 | 0.94
2024-01-23 | 1.61
2024-02-07 | 1.53
2024-02-21 | 1.19
2024-03-06 | 1.17
2024-03-22 | 0.85
2024-04-04 | 1.22
2024-04-16 | 1.10
treatment: MC1
date | Predicted
----------------------
2023-12-27 | 1.43
2024-01-12 | 0.77
2024-01-23 | 0.83
2024-02-07 | 1.27
2024-02-21 | 0.99
2024-03-06 | 0.97
2024-03-22 | 0.71
2024-04-04 | 1.41
2024-04-16 | 0.92
treatment: MC2
date | Predicted
----------------------
2023-12-27 | 1.28
2024-01-12 | 1.00
2024-01-23 | 0.60
2024-02-07 | 1.13
2024-02-21 | 0.88
2024-03-06 | 0.87
2024-03-22 | 0.63
2024-04-04 | 0.93
2024-04-16 | 0.82
treatment: PC0
date | Predicted
----------------------
2023-12-27 | 1.11
2024-01-12 | 0.61
2024-01-23 | 1.04
2024-02-07 | 0.99
2024-02-21 | 0.77
2024-03-06 | 0.75
2024-03-22 | 0.55
2024-04-04 | 0.81
2024-04-16 | 0.82
treatment: PC1
date | Predicted
----------------------
2023-12-27 | 1.69
2024-01-12 | 0.72
2024-01-23 | 3.85
2024-02-07 | 1.50
2024-02-21 | 0.81
2024-03-06 | 1.15
2024-03-22 | 0.84
2024-04-04 | 1.34
2024-04-16 | 1.88
treatment: PC2
date | Predicted
----------------------
2023-12-27 | 0.98
2024-01-12 | 0.70
2024-01-23 | 1.12
2024-02-07 | 0.77
2024-02-21 | 0.68
2024-03-06 | 0.66
2024-03-22 | 0.48
2024-04-04 | 0.82
2024-04-16 | 0.63
Adjusted for:
* site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)
fit_mg <- glmmTMB(
mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
zi = ~ treatment,
family = glmmTMB::lognormal(link = "log"),
data = df_model
)
S(fit_mg) Family: lognormal ( log )
Formula:
mg ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation: ~treatment
Data: df_model
AIC BIC logLik -2*log(L) df.resid
3712.5 3964.0 -1793.3 3586.5 337
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site_block_id (Intercept) 0.08437 0.2905
trt_site_block_id (Intercept) 0.06995 0.2645
Number of obs: 400, groups: site_block_id, 8; trt_site_block_id, 48
Dispersion parameter for lognormal family (): 20.5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.89132 0.19848 19.605 < 2e-16 ***
date2024-01-12 0.19747 0.18513 1.067 0.286123
date2024-01-23 0.31634 0.15943 1.984 0.047241 *
date2024-02-07 0.12808 0.16715 0.766 0.443514
date2024-02-21 0.21117 0.16575 1.274 0.202633
date2024-03-06 0.32937 0.17429 1.890 0.058790 .
date2024-03-22 0.48842 0.16919 2.887 0.003891 **
date2024-04-04 0.55742 0.17043 3.271 0.001073 **
date2024-04-16 0.41760 0.17030 2.452 0.014203 *
treatmentMC1 0.63274 0.21271 2.975 0.002934 **
treatmentMC2 0.25160 0.22258 1.130 0.258313
treatmentPC0 -0.66214 0.25529 -2.594 0.009497 **
treatmentPC1 -0.10683 0.25276 -0.423 0.672552
treatmentPC2 -0.36107 0.28614 -1.262 0.206999
date2024-01-12:treatmentMC1 -0.12862 0.21072 -0.610 0.541599
date2024-01-23:treatmentMC1 -1.38917 0.28870 -4.812 1.5e-06 ***
date2024-02-07:treatmentMC1 -0.50278 0.20987 -2.396 0.016588 *
date2024-02-21:treatmentMC1 -0.54368 0.20679 -2.629 0.008561 **
date2024-03-06:treatmentMC1 -0.64073 0.21213 -3.020 0.002524 **
date2024-03-22:treatmentMC1 -0.64668 0.20587 -3.141 0.001683 **
date2024-04-04:treatmentMC1 -0.66827 0.21603 -3.093 0.001978 **
date2024-04-16:treatmentMC1 -0.78931 0.21531 -3.666 0.000246 ***
date2024-01-12:treatmentMC2 0.10814 0.22548 0.480 0.631515
date2024-01-23:treatmentMC2 0.10015 0.20230 0.495 0.620537
date2024-02-07:treatmentMC2 0.02388 0.21279 0.112 0.910640
date2024-02-21:treatmentMC2 0.10738 0.20869 0.515 0.606875
date2024-03-06:treatmentMC2 0.13378 0.21369 0.626 0.531277
date2024-03-22:treatmentMC2 0.13797 0.20790 0.664 0.506919
date2024-04-04:treatmentMC2 0.05463 0.21106 0.259 0.795775
date2024-04-16:treatmentMC2 -0.30114 0.22976 -1.311 0.189974
date2024-01-12:treatmentPC0 0.25153 0.25958 0.969 0.332538
date2024-01-23:treatmentPC0 0.28432 0.24273 1.171 0.241472
date2024-02-07:treatmentPC0 0.36383 0.25736 1.414 0.157443
date2024-02-21:treatmentPC0 0.41553 0.25202 1.649 0.099185 .
date2024-03-06:treatmentPC0 0.39201 0.25491 1.538 0.124083
date2024-03-22:treatmentPC0 0.28602 0.24940 1.147 0.251441
date2024-04-04:treatmentPC0 0.28631 0.25015 1.145 0.252393
date2024-04-16:treatmentPC0 0.35281 0.26061 1.354 0.175802
date2024-01-12:treatmentPC1 -0.02233 0.26210 -0.085 0.932094
date2024-01-23:treatmentPC1 -0.19236 0.25388 -0.758 0.448642
date2024-02-07:treatmentPC1 -0.22277 0.25231 -0.883 0.377271
date2024-02-21:treatmentPC1 -0.18247 0.25305 -0.721 0.470860
date2024-03-06:treatmentPC1 -0.13884 0.25573 -0.543 0.587192
date2024-03-22:treatmentPC1 -0.16120 0.25385 -0.635 0.525414
date2024-04-04:treatmentPC1 -0.19426 0.26067 -0.745 0.456122
date2024-04-16:treatmentPC1 -0.26135 0.26897 -0.972 0.331209
date2024-01-12:treatmentPC2 0.37069 0.31011 1.195 0.231953
date2024-01-23:treatmentPC2 -0.09502 0.27637 -0.344 0.730988
date2024-02-07:treatmentPC2 -0.21252 0.29873 -0.711 0.476840
date2024-02-21:treatmentPC2 -0.04440 0.29081 -0.153 0.878665
date2024-03-06:treatmentPC2 -0.09098 0.30691 -0.296 0.766902
date2024-03-22:treatmentPC2 -0.05387 0.29106 -0.185 0.853159
date2024-04-04:treatmentPC2 -0.08550 0.29353 -0.291 0.770830
date2024-04-16:treatmentPC2 -0.21145 0.30613 -0.691 0.489735
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.745 2462.968 -0.008 0.994
treatmentMC1 -4.190 19409.569 0.000 1.000
treatmentMC2 -4.202 19239.776 0.000 1.000
treatmentPC0 -4.212 19067.084 0.000 1.000
treatmentPC1 -4.184 19493.173 0.000 1.000
treatmentPC2 15.586 2462.969 0.006 0.995
simulationOutput <- simulateResiduals(fittedModel = fit_mg, n = 500)
plot(simulationOutput)Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
: Fitting terminated with step failure - check results carefully
testUniformity(simulationOutput)
Asymptotic one-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.0445, p-value = 0.4067
alternative hypothesis: two-sided
testOutliers(simulationOutput)
DHARMa outlier test based on exact binomial test with approximate
expectations
data: simulationOutput
outliers at both margin(s) = 1, observations = 400, p-value = 1
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
6.329252e-05 1.384977e-02
sample estimates:
frequency of outliers (expected: 0.00399201596806387 )
0.0025
testDispersion(simulationOutput)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 1.0908, p-value = 0.612
alternative hypothesis: two.sided
testZeroInflation(simulationOutput)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.0941, p-value = 1
alternative hypothesis: two.sided
library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_mg, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") +
theme_bw() +
labs(
x = "Date",
y = "Mg",
title = "Predicted Values of Mg",
color = "Treatment"
)print(pred_effects)# Predicted values of mg
treatment: MC0
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 48.98 | 33.19, 72.26
2024-01-12 | 59.67 | 41.83, 85.12
2024-01-23 | 67.20 | 48.77, 92.58
2024-02-07 | 55.67 | 39.96, 77.55
2024-02-21 | 60.49 | 43.37, 84.36
2024-03-06 | 68.08 | 48.83, 94.92
2024-03-22 | 79.82 | 57.83, 110.17
2024-04-04 | 85.52 | 61.87, 118.22
2024-04-16 | 74.36 | 53.75, 102.88
treatment: MC1
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 92.21 | 66.66, 127.56
2024-01-12 | 98.78 | 72.47, 134.65
2024-01-23 | 31.54 | 18.42, 54.01
2024-02-07 | 63.39 | 45.70, 87.93
2024-02-21 | 66.12 | 48.00, 91.10
2024-03-06 | 67.54 | 49.15, 92.81
2024-03-22 | 78.71 | 57.72, 107.34
2024-04-04 | 82.53 | 59.55, 114.38
2024-04-16 | 63.58 | 45.95, 87.99
treatment: MC2
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 62.99 | 44.42, 89.31
2024-01-12 | 85.50 | 62.61, 116.77
2024-01-23 | 95.53 | 70.42, 129.58
2024-02-07 | 73.32 | 53.58, 100.34
2024-02-21 | 86.61 | 63.86, 117.48
2024-03-06 | 100.09 | 74.39, 134.66
2024-03-22 | 117.84 | 88.03, 157.73
2024-04-04 | 116.16 | 86.53, 155.94
2024-04-16 | 70.77 | 50.38, 99.39
treatment: PC0
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 25.26 | 16.40, 38.90
2024-01-12 | 39.57 | 27.40, 57.16
2024-01-23 | 46.05 | 32.79, 64.68
2024-02-07 | 41.31 | 29.20, 58.45
2024-02-21 | 47.27 | 33.81, 66.09
2024-03-06 | 51.96 | 37.45, 72.11
2024-03-22 | 54.80 | 39.63, 75.76
2024-04-04 | 58.73 | 42.65, 80.86
2024-04-16 | 54.58 | 39.20, 75.99
treatment: PC1
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 44.01 | 28.78, 67.32
2024-01-12 | 52.44 | 36.28, 75.79
2024-01-23 | 49.82 | 34.42, 72.12
2024-02-07 | 40.04 | 28.05, 57.15
2024-02-21 | 45.29 | 32.09, 63.93
2024-03-06 | 53.25 | 38.18, 74.28
2024-03-22 | 61.05 | 44.13, 84.46
2024-04-04 | 63.28 | 45.15, 88.71
2024-04-16 | 51.46 | 35.63, 74.31
treatment: PC2
date | Predicted | 95% CI
--------------------------------------
2023-12-27 | 34.13 | 20.75, 56.15
2024-01-12 | 60.24 | 41.87, 86.68
2024-01-23 | 42.59 | 29.66, 61.16
2024-02-07 | 31.37 | 21.44, 45.89
2024-02-21 | 40.33 | 28.26, 57.55
2024-03-06 | 43.32 | 29.98, 62.60
2024-03-22 | 52.71 | 37.58, 73.94
2024-04-04 | 54.72 | 38.69, 77.38
2024-04-16 | 41.95 | 29.13, 60.41
Adjusted for:
* site_block_id = NA (population-level)
* trt_site_block_id = NA (population-level)
fit_p <- glmmTMB(
p ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id),
zi = ~ treatment,
family = glmmTMB::lognormal(link = "log"),
data = df_model
)
S(fit_p) Family: lognormal ( log )
Formula:
p ~ date * treatment + (1 | site_block_id) + (1 | trt_site_block_id)
Zero inflation: ~treatment
Data: df_model
AIC BIC logLik -2*log(L) df.resid
-807.8 -556.2 466.9 -933.8 338
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site_block_id (Intercept) 0.003353 0.05790
trt_site_block_id (Intercept) 0.007119 0.08438
Number of obs: 401, groups: site_block_id, 8; trt_site_block_id, 48
Dispersion parameter for lognormal family (): 0.0619
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.538415 0.151526 -10.153 < 2e-16 ***
date2024-01-12 -0.880465 0.265132 -3.321 0.000897 ***
date2024-01-23 -0.808236 0.196075 -4.122 3.75e-05 ***
date2024-02-07 -1.320144 0.235561 -5.604 2.09e-08 ***
date2024-02-21 -1.515269 0.306446 -4.945 7.63e-07 ***
date2024-03-06 -1.571143 0.262633 -5.982 2.20e-09 ***
date2024-03-22 -1.364637 0.268976 -5.073 3.91e-07 ***
date2024-04-04 -0.891213 0.251338 -3.546 0.000391 ***
date2024-04-16 -1.407929 0.265990 -5.293 1.20e-07 ***
treatmentMC1 -0.099689 0.195277 -0.511 0.609700
treatmentMC2 -0.158669 0.190661 -0.832 0.405295
treatmentPC0 -0.639774 0.262584 -2.436 0.014832 *
treatmentPC1 0.121799 0.199394 0.611 0.541303
treatmentPC2 -0.009262 0.239372 -0.039 0.969134
date2024-01-12:treatmentMC1 0.409739 0.339480 1.207 0.227447
date2024-01-23:treatmentMC1 -0.091786 0.317385 -0.289 0.772433
date2024-02-07:treatmentMC1 -0.051902 0.348987 -0.149 0.881774
date2024-02-21:treatmentMC1 0.056781 0.392339 0.145 0.884929
date2024-03-06:treatmentMC1 0.161877 0.373965 0.433 0.665111
date2024-03-22:treatmentMC1 0.052483 0.418521 0.125 0.900205
date2024-04-04:treatmentMC1 -0.772945 0.369608 -2.091 0.036505 *
date2024-04-16:treatmentMC1 -0.220208 0.418549 -0.526 0.598804
date2024-01-12:treatmentMC2 0.407385 0.314855 1.294 0.195707
date2024-01-23:treatmentMC2 -0.268135 0.294249 -0.911 0.362162
date2024-02-07:treatmentMC2 -0.014808 0.355481 -0.042 0.966772
date2024-02-21:treatmentMC2 0.077202 0.399456 0.193 0.846749
date2024-03-06:treatmentMC2 0.129368 0.406547 0.318 0.750325
date2024-03-22:treatmentMC2 -0.216069 0.413588 -0.522 0.601374
date2024-04-04:treatmentMC2 -0.234763 0.368105 -0.638 0.523630
date2024-04-16:treatmentMC2 0.299460 0.371436 0.806 0.420116
date2024-01-12:treatmentPC0 0.809856 0.386503 2.095 0.036140 *
date2024-01-23:treatmentPC0 0.070908 0.350091 0.203 0.839494
date2024-02-07:treatmentPC0 0.439496 0.394425 1.114 0.265163
date2024-02-21:treatmentPC0 0.755476 0.442313 1.708 0.087634 .
date2024-03-06:treatmentPC0 0.700833 0.419784 1.670 0.095016 .
date2024-03-22:treatmentPC0 0.378590 0.428899 0.883 0.377397
date2024-04-04:treatmentPC0 0.075187 0.388650 0.193 0.846601
date2024-04-16:treatmentPC0 0.504392 0.389418 1.295 0.195235
date2024-01-12:treatmentPC1 0.277727 0.324591 0.856 0.392207
date2024-01-23:treatmentPC1 -0.516385 0.300661 -1.718 0.085887 .
date2024-02-07:treatmentPC1 -0.113096 0.323280 -0.350 0.726459
date2024-02-21:treatmentPC1 -0.029375 0.421446 -0.070 0.944432
date2024-03-06:treatmentPC1 -0.095028 0.369877 -0.257 0.797243
date2024-03-22:treatmentPC1 -0.285848 0.368584 -0.776 0.438027
date2024-04-04:treatmentPC1 -0.542964 0.349546 -1.553 0.120341
date2024-04-16:treatmentPC1 -0.212736 0.384762 -0.553 0.580330
date2024-01-12:treatmentPC2 0.587871 0.356377 1.650 0.099030 .
date2024-01-23:treatmentPC2 -0.017977 0.315887 -0.057 0.954618
date2024-02-07:treatmentPC2 0.197177 0.348137 0.566 0.571137
date2024-02-21:treatmentPC2 0.012480 0.418305 0.030 0.976199
date2024-03-06:treatmentPC2 0.065331 0.383769 0.170 0.864826
date2024-03-22:treatmentPC2 0.302133 0.397506 0.760 0.447212
date2024-04-04:treatmentPC2 0.319921 0.349066 0.917 0.359403
date2024-04-16:treatmentPC2 0.324149 0.383930 0.844 0.398506
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.53393 0.33244 -4.614 3.95e-06 ***
treatmentMC1 0.51228 0.43136 1.188 0.235
treatmentMC2 0.63784 0.42535 1.500 0.134
treatmentPC0 0.03843 0.45242 0.085 0.932
treatmentPC1 -0.59769 0.51993 -1.150 0.250
treatmentPC2 -0.42968 0.50305 -0.854 0.393
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = fit_p, n = 500)
plot(simulationOutput)testUniformity(simulationOutput)
Asymptotic one-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.087536, p-value = 0.004287
alternative hypothesis: two-sided
testOutliers(simulationOutput)
DHARMa outlier test based on exact binomial test with approximate
expectations
data: simulationOutput
outliers at both margin(s) = 3, observations = 401, p-value = 0.2166
alternative hypothesis: true probability of success is not equal to 0.003992016
95 percent confidence interval:
0.001545489 0.021706793
sample estimates:
frequency of outliers (expected: 0.00399201596806387 )
0.007481297
testDispersion(simulationOutput)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 2.5953, p-value = 0.008
alternative hypothesis: two.sided
testZeroInflation(simulationOutput)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.0005, p-value = 1
alternative hypothesis: two.sided
library(ggeffects)
library(RColorBrewer)
pred_effects <- ggpredict(fit_p, terms = c("date", "treatment"))
plot(pred_effects, colors = "flat") +
theme_bw() +
labs(
x = "Date",
y = "P",
title = "Predicted Values of P",
color = "Treatment"
)